如何对地球表面数据覆盖范围外区域进行球面插值？

问题根源 #

你当前使用的scipy.interpolate.griddata是基于平面笛卡尔坐标系的插值方法，仅能在数据点构成的凸包范围内生成有效结果，超出该范围（比如极地、远离所有数据点的区域）会返回NaN，导致绘图出现空白。地球是球面结构，必须采用基于球面几何的插值方法才能覆盖整个地表。

方案1：使用pyinterp实现专业球面插值 #

pyinterp是专为地理空间数据设计的库，支持高效的球面插值，完美适配全球范围的数据集。

步骤1：安装依赖 #

pip install pyinterp

步骤2：修改后的完整代码 #

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import pyinterp

# 预处理经纬度：将超出[-180, 180]范围的经度修正
lats = np.array([40, 37, 33, -29.7, -34.4, -35.3, 47, -85, 64, 18])
lons = np.array([-105.30, -7.0, -233.7, -54.0, 19.22, 149.0, 47, -85, 64, 18])
lons = np.where(lons > 180, lons - 360, lons)
lons = np.where(lons < -180, lons + 360, lons)
data = np.array([0, 5, 10, 0, 5, 10,7,8,6,8])

# 创建球面插值索引器
interpolator = pyinterp.RTree()
interpolator.packing(
    np.vstack((lons, lats)).T,  # 传入经纬度坐标对（lon, lat）
    data
)

# 生成全球覆盖的网格
xi = np.linspace(-180, 180, 360)
yi = np.linspace(-90, 90, 180)
xi, yi = np.meshgrid(xi, yi)

# 执行球面线性插值，支持全地球范围
zi = interpolator.radial_basis(
    np.vstack((xi.ravel(), yi.ravel())).T,
    method="linear",
    radius=1000000,  # 搜索半径（单位：米），可根据数据密度调整
    num_threads=0
).reshape(xi.shape)

# 绘图逻辑与原代码一致
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=0))
ax.coastlines()
ax.gridlines(draw_labels=True)
ax.set_extent([-180, 180, -90, 90], crs=ccrs.PlateCarree())

# 绘制插值后的颜色图
plt.pcolormesh(xi, yi, zi, shading='auto', cmap='viridis')

# 标记原始数据点
scatter = ax.scatter(lons, lats, c=data, cmap='viridis', s=100,
                     edgecolors='black', linewidths=1,
                     transform=ccrs.PlateCarree())

# 配置颜色条
cbar = plt.colorbar(scatter, label='Number of monkeys (Monkeys /$m^2$)', ax=ax, shrink=0.7, pad=.1)
cbar.mappable.set_clim(0, 10)

plt.title('Spherical Interpolation Colorplot with\nScatter Plot (PlateCarree Projection)')
plt.show()

关键说明 #

• pyinterp的radial_basis方法基于球面距离搜索邻域点，可覆盖包括极地在内的整个地球表面，不会出现空白区域。
• 经度预处理是必要操作，确保所有坐标落在标准的[-180, 180]经度范围内，避免坐标混乱。

方案2：手动转换为三维笛卡尔坐标插值 #

如果不想安装额外库，可以将经纬度转换为三维笛卡尔坐标，在三维空间完成插值后再转回经纬度网格。

代码示例 #

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from scipy.interpolate import griddata

def latlon_to_cartesian(lat, lon):
    """经纬度转三维笛卡尔坐标"""
    lat_rad = np.radians(lat)
    lon_rad = np.radians(lon)
    x = np.cos(lat_rad) * np.cos(lon_rad)
    y = np.cos(lat_rad) * np.sin(lon_rad)
    z = np.sin(lat_rad)
    return x, y, z

def cartesian_to_latlon(x, y, z):
    """三维笛卡尔坐标转经纬度"""
    lat = np.degrees(np.arcsin(z))
    lon = np.degrees(np.arctan2(y, x))
    return lat, lon

# 预处理经度
lats = np.array([40, 37, 33, -29.7, -34.4, -35.3, 47, -85, 64, 18])
lons = np.array([-105.30, -7.0, -233.7, -54.0, 19.22, 149.0, 47, -85, 64, 18])
lons = np.where(lons > 180, lons - 360, lons)
lons = np.where(lons < -180, lons + 360, lons)
data = np.array([0, 5, 10, 0, 5, 10,7,8,6,8])

# 转换原始数据点为三维坐标
x_points, y_points, z_points = latlon_to_cartesian(lats, lons)

# 生成全球网格并转换为三维坐标
xi_lon = np.linspace(-180, 180, 100)
yi_lat = np.linspace(-90, 90, 100)
xi_lon, yi_lat = np.meshgrid(xi_lon, yi_lat)
xi, yi, zi = latlon_to_cartesian(yi_lat, xi_lon)

# 执行三维空间插值
zi_data = griddata((x_points, y_points, z_points), data, (xi, yi, zi), method='linear')

# 绘图逻辑
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=0))
ax.coastlines()
ax.gridlines(draw_labels=True)
ax.set_extent([-180, 180, -90, 90], crs=ccrs.PlateCarree())

plt.pcolormesh(xi_lon, yi_lat, zi_data, shading='auto', cmap='viridis')

scatter = ax.scatter(lons, lats, c=data, cmap='viridis', s=100,
                     edgecolors='black', linewidths=1,
                     transform=ccrs.PlateCarree())

cbar = plt.colorbar(scatter, label='Number of monkeys (Monkeys /$m^2$)', ax=ax, shrink=0.7, pad=.1)
cbar.mappable.set_clim(0, 10)

plt.title('3D Cartesian Interpolation for Spherical Data')
plt.show()

关键说明 #

• 该方法将球面问题转化为三维空间的平面插值，能覆盖全球，但插值精度略低于专业球面插值库，极地附近误差相对明显。

注意事项 #

1. 必须确保所有经度落在[-180, 180]范围内，否则会出现坐标错位问题。
2. 若需要更平滑的插值结果，可尝试pyinterp的gaussian方法，或调整搜索半径参数。

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