如何使用Python与Cartopy绘制按数值填充的美国县域等值区域图

核心思路 #

MetPy内置的USCOUNTIES数据源本身附带每个县域的属性信息（包括FIPS编码、县名、所属州名等），你只需要遍历所有县域几何，匹配你自定义的数值映射，给每个县域分配对应填充色后叠加到图上即可。

具体实现代码 #

以下是在你原有代码基础上修改后的完整可运行版本，包含自定义数值匹配、按值填充的逻辑：

import cartopy
import cartopy.crs as ccrs
import matplotlib as mpl
import matplotlib.pyplot as plt
from metpy.plots import USCOUNTIES
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
import numpy as np

plot_type = 'png'

# ---------------------- 自定义配置部分 ----------------------
# 1. 定义县域到数值的映射，推荐用5位FIPS码做键，避免同名县冲突
# 你也可以替换为"县名,州名"格式的字符串做键，比如"洛杉矶县,加利福尼亚州"
county_value_map = {
    '06037': 2,   # 加利福尼亚州洛杉矶县：居住过
    '17031': 1,   # 伊利诺伊州库克县：到访过
    '48201': 1,   # 得克萨斯州哈里斯县：到访过
    '36061': 2    # 纽约州纽约县：居住过
}

# 2. 定义颜色映射规则，可根据你的数值类型调整
# 示例：0=无记录（透明）、1=到访（浅绿）、2=居住（深绿）
value_color_map = mpl.colors.ListedColormap(['none', '#90EE90', '#006400'])
value_norm = mpl.colors.Normalize(vmin=0, vmax=2)
# 如果是连续数值（比如龙卷风评级、人口），可以用连续配色：
# value_color_map = plt.get_cmap('viridis')
# value_norm = mpl.colors.Normalize(vmin=0, vmax=5) # 对应EF0-EF5龙卷风评级
# -----------------------------------------------------------

borders  = cartopy.feature.BORDERS
states   = cartopy.feature.NaturalEarthFeature(category='cultural', scale='10m', facecolor='none', name='admin_1_states_provinces_lakes')
oceans   = cartopy.feature.OCEAN
lakes    = cartopy.feature.LAKES

mpl.rcParams['figure.figsize'] = (12,10)
water_color = 'lightblue'

fig = plt.figure()
ax = plt.axes(projection=ccrs.LambertConformal(central_longitude=-97.5, central_latitude=38.5, standard_parallels=(38.5,38.5)))
ax.set_extent([-120, -74, 23, 50], ccrs.Geodetic())
ax.coastlines()
ax.add_feature(borders, linestyle='-')
ax.add_feature(states, linewidth=0.50, edgecolor='black')
ax.add_feature(oceans, facecolor=water_color)
ax.add_feature(lakes, facecolor=water_color, linewidth=0.50, edgecolor='black')

# ---------------------- 县域填充逻辑 ----------------------
county_reader = USCOUNTIES.with_scale('500k')
patches = []
fill_colors = []

for county in county_reader.records():
    # 提取县域属性，可根据你的匹配规则调整
    county_key = county.attributes['GEOID'] # 用FIPS码做匹配键
    # 如果要用县名+州名匹配，替换为下面的代码：
    # county_key = f"{county.attributes['NAME']}, {county.attributes['STATE_NAME']}"
    
    # 匹配数值，无匹配默认返回0（透明）
    current_value = county_value_map.get(county_key, 0)
    # 获取县域几何，处理MultiPolygon类型的几何
    geom = county.geometry
    polys = [geom] if geom.geom_type == 'Polygon' else list(geom.geoms)
    
    for poly in polys:
        # 转换坐标到当前投影
        coord_arr = np.array(poly.exterior.coords)
        proj_coords = ccrs.PlateCarree().transform_points(ccrs.PlateCarree(), coord_arr[:, 0], coord_arr[:, 1])
        # 生成多边形patch
        patches.append(Polygon(proj_coords[:, :2], closed=True))
        fill_colors.append(value_color_map(value_norm(current_value)))

# 把所有填充多边形添加到图上，注意要在加县域边界之前添加，避免边界被填充色覆盖
fill_collection = PatchCollection(patches, match_original=True)
fill_collection.set_facecolor(fill_colors)
fill_collection.set_edgecolor('none')
ax.add_collection(fill_collection)
# -----------------------------------------------------------

# 最后加县域边界，保证边界显示在填充色上层
ax.add_feature(USCOUNTIES.with_scale('500k'), linewidth=0.10, edgecolor='black')

# 如果是连续数值可以加 colorbar
# sm = plt.cm.ScalarMappable(cmap=value_color_map, norm=value_norm)
# sm.set_array([])
# plt.colorbar(sm, ax=ax, shrink=0.6, label='你的数值标签')

plt.savefig('./county_map.'+plot_type, dpi=300, bbox_inches='tight')
plt.close()

注意事项 #

• 优先使用5位FIPS编码作为匹配键，美国每个县域有唯一的FIPS编码，完全避免同名县匹配错误的问题
• 如果需要添加阿拉斯加、夏威夷的插图，只需新建两个小尺寸子图，调整子图的set_extent到对应区域，重复上述填充逻辑即可
• 数值为离散分类还是连续值，只需要调整颜色映射（Colormap）和归一化规则（Normalize）即可适配

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