Geopandas轻松处理shp文件的投影转换

        在绘制一张降水的地理分布图时，始终不能将地图边界在图上显示出来，经过一番折腾，最后发现问题出在shp文件的地图投影：该投影下的边界并不是以经纬度（如 [121.22269, 40.94580]）给出，而是以X、Y坐标值给出（如[ 1338807.250, 4516312.500,]），不能显示地图的原因是错将后者作为经纬度。

        查阅了几种改变投影的工具，pyshp，pyproj的Transformer，Proj，CRS等，看的晕头转向，却仍然没有将X、Y坐标正确转换为经纬度坐标。最后是这篇文章（Map projections — Intro to Python GIS documentation）帮助了我，轻松用Geopandas搞定。以下是大体过程：

import geopandas as gpd

# 读入原shp文件
shp = gpd.read_file('../shp_files/NINE_BASINS/liuyu.shp') 

# 利用.crs 属性获取投影信息（实际我们并不需要知道这些详细信息，to_crs会自动使用）
print(shp.crs)  

# 利用.copy()方法复制为新的shp对象，用来变换投影
shp_new = shp.copy()

# 利用.to_crs()方法将原投影转换为epsg关键字所指定的投影，这里4326对应WGS84投影。
shp_new = shp_new.to_crs(epsg=4326) 

# 利用字段值索引其中一个子区域用于绘图（shp文件含9大流域，这里只需要珠江流域）。
shp_zhujiang = shp_new[shp_new['W1102WB0_1']=="珠江流域片"]['geometry']

#----
# 此处略去绘图过程
#----

# 利用add_geometries在等值线图上叠加地图边界shp_zhujiang。
ax.add_geometries(shp_zhujiang, crs=projection, facecolor='none', edgecolor='k', linewidth=2,zorder=1)

绘图结果如下：

注：

1. 所用的9大流域片边界数据来源：中国科学院资源环境科学数据中心。采用的投影方式是："Albers_Conic_Equal_Area"。

2. epsg=4326中的数字是Well Known ID (WKID)，其与投影的对应关系，可以查看地理坐标（Coordinate Systems）和投影坐标（Coordinate Systems）。