最近在写毕业论文,上手练习一下cartopy包。
问题:
正常运行ax.add_geometries代码,但是shp文件在结果图中一直不显示。一开始还以为是投影有问题,代码中设置的是简易圆柱(Plate Carrée) 投影,但后面也尝试了其他投影,依然显示不了。
ax.add_geometries(Reader(shp_path).geometries(), ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=0.65)
运行结果如下:
问题发现:
原来是子图中的绘图范围与Shp的数值范围不匹配!
即:我之前提供的shp文件是投影过的,检查代码如下:
import geopandas as gpd
# 读取Shapefile
gdf = gpd.read_file(shp_path)
# 获取边界
bounds = gdf.total_bounds
print("Shapefile bounds:", bounds)
#结果: [11448677.0 3047834.8 11625463.8 3218060.9 ]
然而,我的绘图范围确是用经纬度表示的,检查代码如下:
xlim = ax.get_xlim()#这个ax是我子图的轴
ylim = ax.get_ylim()
print("Current map view limits:", xlim, ylim)
#结果:Current map view limits: (102.5, 105.0) (27.0, 29.5)
解决:
将shp文件投影更替为地理坐标系即可,我是用的Arcgis的project工具。然后,再检查shp文件的数值范围时,其范围就落入绘图范围中啦!
#同样的代码
import geopandas as gpd
# 读取Shapefile
gdf = gpd.read_file(shp_path)
# 获取边界
bounds = gdf.total_bounds
print("Shapefile bounds:", bounds)
#结果更新: [102.9 27.5 104.7 28.9]
绘图结果:
完结撒花!