Geopandas 的一些常规操作

geopandas添加底图

import contextily as ctx # 这个库是用来添加底图的

fig, ax = plt.subplots(figsize=(10, 15))
manhattan = manhattan.to_crs('EPSG:3857')
ax = manhattan.plot(ax=ax, alpha=0.4, edgecolor='k')
ax.axis('off')
centroid = manhattan.geometry.centroid
for row in manhattan.itertuples():
    ax.text(centroid.x.tolist()[row.Index], centroid.y.tolist()[row.Index], row.Index, ha='center', va='center')

# 添加底图、不过加上后，这个图片显示的有点慢，可以更换url，以更换不同风格的底图
ctx.add_basemap(ax, 
				source='https://d.basemaps.cartocdn.com/light_nolabels/{z}/{x}/{y}.png',
                zoom=13)
plt.savefig('images/地理区域图_排序id.png', dpi=300)

显示如图：

显示每个区域的标签

fig, ax = plt.subplots(1, figsize = (10,10))
cd.plot(ax=ax, figsize=(10, 10), linewidth=0.8, edgecolor='gray', legend = True, color='y')
cd.apply(lambda x:ax.annotate(text=x.ID, xy=x.geometry.centroid.coords[0], ha='center', color='black'), axis=1)
ax.axis('off')

效果图，但它好像会打印一些不需要的东西（如下面红框中的内容），还不知道怎么取消

两个GeoDataFram的连接

manhattan = pd.merge(ny, region, on=['location_id']) #其实就是pandas的合并操作

将DataFrame转换成GeoDataFrame

首先需要将经纬度转换成geometry形式的点、线、面

from shapely import geometry

data['PULocation'] = data.apply(lambda row:geometry.Point(row['SLon'], row['SLat']), axis=1)
# 下面这两步好像有点多余
point1 = gpd.GeoSeries(data['PULocation'], crs='EPSG:4326') 
data['PULocation'] = point1
pu = gpd.GeoDataFrame(pu, geometry='PULocation', crs='EPSG:4326')

做空间连接

如选出data中的点在cd面中的哪个区域

# op参数可以选择
join = gpd.sjoin(data, cd, how='left', op='within')

op：用于设定拓扑判断的规则，'intersects’代表相交，即几何对象之间存在共有的边或内部点；'contains’代表包含，即一个几何对象至少有一个点位于另一个几何对象内部，且其本身没有任何点落在另一个结几何对象的外部；‘within’表示在内部，是’contains’的相反情况，即左表被右表矢量’contains’

地理坐标与平面坐标转换

即经纬度坐标 与 xy坐标的转换
本次使用的是osmnx这个库

road = gpd.read_file("D:/Code/data_NYC/shapefile/shape/roads.shp")
road.crs = "EPSG:4326"  # 给定指定一个坐标系不然会报错，这个坐标系我是用QGIS打开后看到的
b = ox.project_gdf(road) 

这样就可以将原来的经纬度转换成xy坐标，转换后如图


xy坐标转换成经纬度：

# b.crs = "EPSG:32618"  # 有时候可能需要也给它指定一个坐标系
ox.project_gdf(b, to_latlong=True) # to_latlong参数默认是False，所以此处需要指定一下

然后就转换回来了。
再用QGIS打开，单位就从度变成了米，长度也正常了