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打开,单位就从度变成了米,长度也正常了