如何用Geopandas为绘图添加美国州界？求改脚本及Shapefile推荐

问题描述 #

需要修改Python脚本，通过Geopandas为Basemap绘制的降水异常图添加美国州界，要求不能使用map.drawstates()方法（该方法会绘制其他国家的州），同时推荐合适的美国州界Shapefile。当前绘图已满足预期，仅缺少美国州界。

关键修改点 #

你的代码无法显示州界的核心原因是：Basemap使用的是投影后的平面坐标，而Shapefile默认是WGS84经纬度坐标，两者坐标系不匹配。需要先将Shapefile的几何数据转换为Basemap的投影坐标系，同时过滤掉非本土州避免干扰。

具体修改步骤 #

• 转换Shapefile坐标系：提取Basemap的投影参数，用to_crs()方法将Geopandas的GeoDataFrame转换为对应投影
• 过滤非本土州：通过州FIPS代码或名称筛选美国本土48州+DC，排除阿拉斯加、夏威夷、海外领地
• 正确指定绘图轴：Basemap绘图时要指定ax=ax（Matplotlib的轴对象），而非直接传入Basemap实例

推荐的美国州界Shapefile #

• US Census Bureau 州界数据：官方权威数据，包含完整的边界属性，推荐使用最新版的tl_20XX_us_state系列，该数据包含所有州及领地，可通过属性筛选本土区域。
• Natural Earth 美国州界：轻量化数据，适合快速绘图，可直接下载美国本土州界的Shapefile，文件体积小，加载速度快。

修改后的完整代码 #

from netCDF4 import Dataset as NetCDFFile
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.colors as mcolors
import geopandas as gpd

# 读取美国州界Shapefile
us_states = gpd.read_file('C:/Users/scapa/Downloads/tl_2012_us_state/tl_2012_us_state.shp')

# 过滤美国本土48州+DC（排除阿拉斯加、夏威夷、波多黎各等）
# 用FIPS代码筛选：阿拉斯加02，夏威夷15，波多黎各72，其余本土州保留
exclude_fips = ['02', '15', '72']
us_mainland = us_states[~us_states['STATEFP'].isin(exclude_fips)]

nc = NetCDFFile('C:/Users/scapa/Downloads/nclwAXbfFo6UE.nc')
lat = nc.variables['lat'][:]
lon = nc.variables['lon'][:]
prate = nc.variables['VAR'][:]

fig, ax = plt.subplots(figsize=(10, 8))
# 初始化Basemap，指定ax参数绑定Matplotlib轴
map = Basemap(llcrnrlon=230., llcrnrlat=10., urcrnrlon=305., urcrnrlat=55., ax=ax)
lons, lats = np.meshgrid(lon, lat)
x, y = map(lons, lats)

# 设置色阶和自定义配色
levels = np.arange(-3.6, 3.6, 0.4)
cmap = plt.get_cmap('BrBG', len(levels) - 1)
cmaplist = [(1, 1, 1, 1) if -0.4 <= val <= 0.4 else cmap((val - levels.min()) / (levels.max() - levels.min())) for val in levels]
custom_cmap = mcolors.LinearSegmentedColormap.from_list('custom', cmaplist, len(levels))

# 将Shapefile转换为Basemap的投影坐标系
# 提取Basemap的proj4字符串
map_proj = map.proj4string
us_mainland_proj = us_mainland.to_crs(map_proj)

# 在Matplotlib轴上绘制州界
us_mainland_proj.boundary.plot(ax=ax, linewidth=0.5, color='black')

# 绘制降水异常填色图
cs = map.contourf(x, y, prate, levels, cmap=custom_cmap)
cbar = plt.colorbar(cs, orientation='horizontal', ticks=levels, pad=0.07)
cbar.set_label('Surface Precipitation Rate Anomalies (mm/day)')
cbar.ax.tick_params(labelsize=8)

# 绘制其他地图元素
map.drawcoastlines()
map.drawparallels(np.arange(10, 60, 10), labels=[1, 1, 1, 1], linewidth=0.5)
map.drawmeridians(np.arange(240, 310, 20), labels=[1, 1, 0, 1], linewidth=0.5)

# 添加标题和副标题
plt.title('El Niño Winter Precipitation Anomalies', fontsize=15, y=1.05)
plt.text(0.5, 1.02, 'Winters of 1957-58, 1965-66, 1972-73, 1982-83, 1987-88, 1991-92, 1997-98, 2009-10, 2015-16', 
         horizontalalignment='center', fontsize=9, fontstyle='italic', transform=ax.transAxes)

plt.savefig('era5.precip.anom.states.png', dpi=600)
plt.show()

说明 #

• 代码中通过STATEFP（州FIPS代码）筛选本土州，确保只显示美国大陆区域的州界
• 转换坐标系后，Shapefile的几何数据才能和Basemap的投影坐标匹配，正确显示在图上
• 绘制州界时指定color='black'可以让边界更清晰，可根据需求调整颜色和线宽

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