geopandas：数据读取、坐标设置与转换、空间查询连接、数据输出

导入相关包

import geopandas as gp
import pandas as pd

数据读取

polygonShp=gp.read_file("./shp/globeland30-10.shp")
polygonShp

	OBJECTID	Id	gridcode	Shape_Leng	Shape_Area	geometry
0	1	1	10	1740.0	59400.000000	POLYGON ((321652.212 3472141.992, 321562.212 3...
1	2	2	10	1020.0	27900.000001	POLYGON ((320602.212 3472111.992, 320572.212 3...
2	3	3	10	1080.0	20700.000000	POLYGON ((320842.212 3472111.992, 320812.212 3...
3	4	4	10	120.0	900.000000	POLYGON ((318442.212 3472021.992, 318412.212 3...
4	5	5	10	120.0	900.000000	POLYGON ((318502.212 3471991.992, 318472.212 3...
...	...	...	...	...	...	...
2998	2999	2999	10	2040.0	107100.000000	POLYGON ((341932.212 3376801.992, 341752.212 3...
2999	3000	3000	10	5100.0	202500.000000	POLYGON ((342412.212 3376171.992, 342382.212 3...
3000	3001	3001	10	5280.0	337500.000001	POLYGON ((342052.212 3375511.992, 341932.212 3...
3001	3002	3002	10	1380.0	60300.000001	POLYGON ((340822.212 3375061.992, 340552.212 3...
3002	3003	3003	10	1620.0	74700.000000	POLYGON ((341812.212 3374521.992, 341692.212 3...

3003 rows × 6 columns

基于带有坐标的Pandas DataFrame创建GeoDataFrame

pointcsv=pd.read_csv("./datadata-Test.csv")
pointcsv

	OBJECTID	Field1	Field2	Field3
0	1	109.353080	30.492034	-9.465
1	2	109.346870	30.491565	3.000
2	3	109.346467	30.492398	8.000
3	4	109.346567	30.492413	4.000
4	5	109.351589	30.494359	39.000
...	...	...	...	...
8429627	8429628	109.047030	31.368869	0.000
8429628	8429629	109.047196	31.368895	0.000
8429629	8429630	109.046931	31.368981	0.000
8429630	8429631	109.047052	31.369000	0.000
8429631	8429632	109.046820	31.369092	0.000

8429632 rows × 4 columns

pointShp = gp.GeoDataFrame(pointcsv, geometry=gp.points_from_xy(pointcsv.Field1,pointcsv.Field2))
pointShp

	OBJECTID	Field1	Field2	Field3	geometry
0	1	109.353080	30.492034	-9.465	POINT (109.35308 30.49203)
1	2	109.346870	30.491565	3.000	POINT (109.34687 30.49157)
2	3	109.346467	30.492398	8.000	POINT (109.34647 30.49240)
3	4	109.346567	30.492413	4.000	POINT (109.34657 30.49241)
4	5	109.351589	30.494359	39.000	POINT (109.35159 30.49436)
...	...	...	...	...	...
8429627	8429628	109.047030	31.368869	0.000	POINT (109.04703 31.36887)
8429628	8429629	109.047196	31.368895	0.000	POINT (109.04720 31.36889)
8429629	8429630	109.046931	31.368981	0.000	POINT (109.04693 31.36898)
8429630	8429631	109.047052	31.369000	0.000	POINT (109.04705 31.36900)
8429631	8429632	109.046820	31.369092	0.000	POINT (109.04682 31.36909)

8429632 rows × 5 columns

polygonShp.plot()

GeoDataFrame设置坐标系与进行坐标转换

# 设置坐标系
pointShp.set_crs(epsg=4326,inplace=True)
pointShp.to_crs(epsg=32649,inplace=True)

两个GeoDataFrame间进行空间连接

join_inner_df = gp.sjoin(pointShp, polygonShp, how='inner', op='intersects')
join_inner_df

	OBJECTID_left	Field1	Field2	Field3	geometry	index_right	OBJECTID_right	Id	gridcode	Shape_Leng	Shape_Area
5	6	109.352016	30.494426	39.0	POINT (341838.913 3374728.302)	3002	3003	3003	10	1620.0	74700.000000
345	346	109.340178	30.497293	1.0	POINT (340707.335 3375062.731)	3001	3002	3002	10	1380.0	60300.000001
346	347	109.340282	30.497310	7.0	POINT (340717.363 3375064.410)	3001	3002	3002	10	1380.0	60300.000001
347	348	109.340370	30.497324	6.0	POINT (340725.746 3375065.815)	3001	3002	3002	10	1380.0	60300.000001
348	349	109.340463	30.497338	5.0	POINT (340734.720 3375067.319)	3001	3002	3002	10	1380.0	60300.000001
...	...	...	...	...	...	...	...	...	...	...	...
8428437	8428438	109.080459	31.370270	0.0	POINT (317445.086 3472230.637)	6	7	7	10	1860.0	108900.000000
8428553	8428554	109.079964	31.370321	0.0	POINT (317398.122 3472237.057)	6	7	7	10	1860.0	108900.000000
8428554	8428555	109.080104	31.370343	-1.0	POINT (317411.484 3472239.253)	6	7	7	10	1860.0	108900.000000
8428671	8428672	109.079921	31.370442	-1.0	POINT (317394.240 3472250.568)	6	7	7	10	1860.0	108900.000000
8428672	8428673	109.080068	31.370465	-1.0	POINT (317408.250 3472252.871)	6	7	7	10	1860.0	108900.000000

4393811 rows × 11 columns

结果输出

join_inner_df.to_csv("./result.csv",columns=['OBJECTID_left','Field1','Field2','Field3'],index=False)