同属EPSG:26917的shp与DEM坐标不匹配，如何统一坐标系使用zonal_stats？

问题根源 #

你的矢量shp文件存在CRS元数据标注错误：
EPSG:26917是UTM 17N投影坐标系，单位为米，正常x轴范围为十万级数值，和你DEM的坐标范围一致。但你shp的x轴范围为-81.9~-81.5，明显是地理坐标系（单位为十进制度）的经度数值，说明shp实际存储的是NAD83地理坐标系（EPSG:4269）的坐标，只是被错误标注成了EPSG:26917。二者坐标范围完全不重叠，所以rasterstats.zonal_stats返回空值。

解决步骤 #

• 第一步：修正shp的CRS元数据，先给矢量数据设置正确的原生坐标系，不要直接重投影

# 覆盖错误的CRS标注，设置为实际对应的NAD83地理坐标系
gdf = gdf.set_crs(epsg=4269, allow_override=True)

• 第二步：将修正后的矢量数据重投影到和DEM一致的EPSG:26917投影坐标系

gdf_proj = gdf.to_crs(epsg=26917)

• 第三步：验证坐标范围匹配
  打印gdf_proj.total_bounds，查看x轴范围是否和DEM的420000~500000一致，确认匹配后即可正常同图绘制，调用rasterstats.zonal_stats也会返回正确统计结果。

内容的提问来源于stack exchange，提问作者Xudong