如何从全球NetCDF气温异常数据中提取单个国家的对应数据

你已经成功读取NASA的全球气温异常NetCDF数据并绘制了全球分布图，要提取德国、法国、尼泊尔、印度等单个国家的数据并单独绘图，可以按以下步骤实现：

步骤1：准备依赖库 #

首先安装需要的地理数据处理库：

pip install geopandas cartopy

步骤2：加载国家边界数据 #

使用geopandas加载内置的自然地球低分辨率国家边界数据集：

import geopandas as gpd
import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset
import cartopy.crs as ccrs

# 读取国家边界数据
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

步骤3：筛选目标国家并提取边界 #

从数据集中筛选出你需要的国家：

# 定义目标国家列表，名称需与数据集内的名称匹配
target_countries = ['Germany', 'France', 'Nepal', 'India']

# 筛选目标国家的边界
country_shapes = world[world['name'].isin(target_countries)]

步骤4：提取对应国家的气温异常数据 #

将气温数据网格点与国家边界匹配，提取落在国家范围内的数据：

# 读取NetCDF数据（复用你已有的代码）
data = Dataset("../data/amaps_robinson_1000km.nc")
lons = data.variables["lon"][:]
lats = data.variables["lat"][:]
temp_anomaly = data.variables["TEMPANOMALY"][:]

# 创建经纬度网格的二维数组
lon_grid, lat_grid = np.meshgrid(lons, lats)

# 将网格点转换为geopandas的点对象
from shapely.geometry import Point
points = [Point(lon, lat) for lon, lat in zip(lon_grid.flatten(), lat_grid.flatten())]
gdf_points = gpd.GeoDataFrame({'geometry': points}, crs="EPSG:4326")

以德国为例，提取单个国家的数据：

# 选择目标国家
country = country_shapes[country_shapes['name'] == 'Germany'].iloc[0]

# 判断每个网格点是否在国家边界内
mask = gdf_points.within(country['geometry'])
# 将掩码转换为和气温数据相同的形状
mask = mask.values.reshape(temp_anomaly.shape)

# 提取国家范围内的气温异常数据，超出范围设为NaN
country_temp = np.where(mask, temp_anomaly, np.nan)

步骤5：绘制单个国家的气温异常图 #

使用cartopy绘制只包含目标国家的分布图：

fig = plt.figure(figsize=(8, 6))
ax = plt.axes(projection=ccrs.PlateCarree())

# 设置地图范围为目标国家的边界范围
bbox = country['geometry'].bounds
ax.set_extent([bbox[0], bbox[2], bbox[1], bbox[3]], crs=ccrs.PlateCarree())

# 绘制国家边界
ax.add_geometries([country['geometry']], crs=ccrs.PlateCarree(), edgecolor='black', facecolor='none', linewidth=1.5)

# 绘制气温异常填充图
clevs = np.arange(-4, 5)
cmap = "coolwarm"
plt.contourf(lons, lats, country_temp, clevs, transform=ccrs.PlateCarree(), cmap=cmap)

# 添加标题和颜色条
plt.title(f"Germany Temperature Anomaly (°C) 2021 vs 1951-1980")
cb = plt.colorbar(orientation="horizontal", pad=0.05, shrink=0.8)
cb.set_label("°C", size=12, rotation=0, labelpad=10)

plt.show()

批量处理多个国家 #

如果要批量处理所有目标国家，可通过循环实现：

clevs = np.arange(-4, 5)
cmap = "coolwarm"

for country_name in target_countries:
    # 筛选国家
    country = country_shapes[country_shapes['name'] == country_name].iloc[0]
    # 创建掩码并提取数据
    mask = gdf_points.within(country['geometry']).values.reshape(temp_anomaly.shape)
    country_temp = np.where(mask, temp_anomaly, np.nan)
    
    # 绘图
    fig = plt.figure(figsize=(8, 6))
    ax = plt.axes(projection=ccrs.PlateCarree())
    bbox = country['geometry'].bounds
    ax.set_extent([bbox[0], bbox[2], bbox[1], bbox[3]], crs=ccrs.PlateCarree())
    ax.add_geometries([country['geometry']], crs=ccrs.PlateCarree(), edgecolor='black', facecolor='none', linewidth=1.5)
    
    plt.contourf(lons, lats, country_temp, clevs, transform=ccrs.PlateCarree(), cmap=cmap)
    plt.title(f"{country_name} Temperature Anomaly (°C) 2021 vs 1951-1980")
    cb = plt.colorbar(orientation="horizontal", pad=0.05, shrink=0.8)
    cb.set_label("°C", size=12, rotation=0, labelpad=10)
    
    plt.savefig(f"{country_name}_temp_anomaly.png", dpi=150)
    plt.close()

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