Pykrige结合Basemap地图插值优化与功能问题排查

原有核心功能代码 #

数据处理与网格生成函数 #

def get_data(df):
    return {
        "lons": df['Longitude'],
        "lats": df['Latitude'],
        "values": df['O18'],
        "alts": df['Altitude'],
    }

def extend_data(data):
    return {
        "lons": np.concatenate([np.array([lon-360 for lon in data["lons"]]), data["lons"], np.array([lon+360 for lon in data["lons"]])]),
        "lats":  np.concatenate([data["lats"], data["lats"], data["lats"]]),
        "values":  np.concatenate([data["values"], data["values"], data["values"]]),
        "alts": np.concatenate([data["alts"], data["alts"], data["alts"]]),
    }

def generate_grid(data, basemap, delta=1):
    grid = {
        'lon': np.arange(-180, 180, delta),
        'lat': np.arange(np.amin(data["lats"]), np.amax(data["lats"]), delta) # 不对极点做外推
    }
    grid["x"], grid["y"] = np.meshgrid(grid["lon"], grid["lat"])
    grid["x"], grid["y"] = basemap(grid["x"], grid["y"])
    return grid

插值函数 #

def interpolate(data, grid):
    UK = UniversalKriging(
        data["lons"],
        data["lats"],
        data["values"],
        variogram_model='exponential',
        specified_drift = data["alts"], 
    )
    return UK.execute("grid", grid["lon"], grid["lat"])

绘图函数 #

def plot_mesh_data(interpolation, grid, basemap):
    colormesh = basemap.contourf(grid["x"], grid["y"],  interpolation,100, cmap='jet', )
    color_bar = basemap.colorbar(colormesh,location='bottom',pad="10%") 

分问题修复方案 #

1. 欧洲区域插值精度不足问题 #

问题原因

• 网格分辨率过低：默认1度间隔的网格无法匹配欧洲区域密集的监测站点分布，细节被过度平滑
• 插值模型参数未优化：默认变差函数滞后距数量不足，全局使用同一套参数未考虑区域站点密度差异
• 渲染层级不足：等值面分段数过少导致边界过渡粗糙

修复代码

# 1. 调整网格分辨率为0.5度，欧洲局部区域可单独裁剪设置0.25度进一步提升精度
def generate_grid(data, basemap, delta=0.5):
    grid = {
        'lon': np.arange(-180, 180, delta),
        'lat': np.arange(np.amin(data["lats"]), np.amax(data["lats"]), delta)
    }
    grid["x"], grid["y"] = np.meshgrid(grid["lon"], grid["lat"])
    grid["x"], grid["y"] = basemap(grid["x"], grid["y"])
    return grid

# 2. 优化克里金插值参数
from pykrige.uk import UniversalKriging
def interpolate(data, grid):
    UK = UniversalKriging(
        data["lons"],
        data["lats"],
        data["values"],
        variogram_model='exponential',
        specified_drift = data["alts"], # 保留高程作为漂移项
        nlags=60, # 提升滞后距数量匹配站点密度
        variogram_parameters={'sill': 0.9, 'range': 15, 'nugget': 0.1}, # 可根据半方差图拟合结果调整
        pseudo_inv=True
    )
    return UK.execute("grid", grid["lon"], grid["lat"])

# 3. 提升等值面渲染精度
def plot_mesh_data(interpolation, grid, basemap):
    colormesh = basemap.contourf(
        grid["x"], grid["y"], interpolation,
        levels=120, cmap='RdBu_r', antialiased=True
    )
    color_bar = basemap.colorbar(colormesh, location='bottom', pad="10%")

2. 南极区域插值权重提升问题 #

实现逻辑

PyKrige不直接支持样本权重参数，可通过对南极区域（纬度<-60°）的观测点做重复采样，人为提升该区域点在插值计算中的影响权重，权重倍数可根据展示效果调整（示例设置为4倍权重），在数据扩展步骤前加入权重处理即可。

修复代码

def weight_antarctica_data(data, weight=4, lat_threshold=-60):
    # 筛选南极区域点
    ant_idx = np.where(data["lats"] < lat_threshold)[0]
    # 按权重重复采样
    ant_lons = np.tile(data["lons"][ant_idx], weight-1)
    ant_lats = np.tile(data["lats"][ant_idx], weight-1)
    ant_values = np.tile(data["values"][ant_idx], weight-1)
    ant_alts = np.tile(data["alts"][ant_idx], weight-1) if "alts" in data else np.array([])
    # 合并回原数据
    data["lons"] = np.concatenate([data["lons"], ant_lons])
    data["lats"] = np.concatenate([data["lats"], ant_lats])
    data["values"] = np.concatenate([data["values"], ant_values])
    if "alts" in data:
        data["alts"] = np.concatenate([data["alts"], ant_alts])
    return data

调用时在extend_data前执行base_data = weight_antarctica_data(base_data)即可生效。 #

3. 格陵兰区域掩膜报错问题 #

问题原因

• 拼写错误：readshapefile注册的属性名是greenland，代码中误写为greendland
• 逻辑错误：格陵兰行政区划shp文件无nombre字段，不存在Selva分类判断，多余判断导致无法生成多边形
• 依赖缺失：未导入Polygon和PatchCollection类

修复代码

首先补充导入依赖：

from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection

修正掩膜函数：

def mask_greenland(axes, basemap, shp_path):
    basemap.readshapefile(shp_path, 'greenland', drawbounds=False)
    patches = []
    # 直接遍历所有格陵兰边界形状，不需要额外属性判断
    for shape in basemap.greenland:
        patches.append(Polygon(np.array(shape), True))
    # 生成掩膜覆盖层，zorder设置高于插值层即可实现遮挡
    mask = PatchCollection(patches, facecolor='white', edgecolor='k', linewidths=0.5, zorder=3)
    axes.add_collection(mask)

调用时传入shp文件的实际路径即可。 #

完整修正后运行代码 #

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
from pykrige.uk import UniversalKriging

def load_data():
    df = pd.read_csv(r"你的数据文件路径.csv")
    return df

def get_data(df):
    return {
        "lons": df['Longitude'].values,
        "lats": df['Latitude'].values,
        "values": df['O18a'].values,
        "alts": df['Altitude'].values if 'Altitude' in df.columns else np.zeros(len(df))
    }

def weight_antarctica_data(data, weight=4, lat_threshold=-60):
    ant_idx = np.where(data["lats"] < lat_threshold)[0]
    ant_lons = np.tile(data["lons"][ant_idx], weight-1)
    ant_lats = np.tile(data["lats"][ant_idx], weight-1)
    ant_values = np.tile(data["values"][ant_idx], weight-1)
    ant_alts = np.tile(data["alts"][ant_idx], weight-1)
    data["lons"] = np.concatenate([data["lons"], ant_lons])
    data["lats"] = np.concatenate([data["lats"], ant_lats])
    data["values"] = np.concatenate([data["values"], ant_values])
    data["alts"] = np.concatenate([data["alts"], ant_alts])
    return data

def extend_data(data):
    return {
        "lons": np.concatenate([np.array([lon-360 for lon in data["lons"]]), data["lons"], np.array([lon+360 for lon in data["lons"]])]),
        "lats":  np.concatenate([data["lats"], data["lats"], data["lats"]]),
        "values":  np.concatenate([data["values"], data["values"], data["values"]]),
        "alts": np.concatenate([data["alts"], data["alts"], data["alts"]]),
    }

def generate_grid(data, basemap, delta=0.5):
    grid = {
        'lon': np.arange(-180, 180, delta),
        'lat': np.arange(np.amin(data["lats"]), np.amax(data["lats"]), delta)
    }
    grid["x"], grid["y"] = np.meshgrid(grid["lon"], grid["lat"])
    grid["x"], grid["y"] = basemap(grid["x"], grid["y"])
    return grid

def interpolate(data, grid):
    UK = UniversalKriging(
        data["lons"],
        data["lats"],
        data["values"],
        variogram_model='exponential',
        specified_drift = data["alts"],
        nlags=60,
        pseudo_inv=True,
        verbose=True,
    )
    return UK.execute("grid", grid["lon"], grid["lat"])

def prepare_map_plot():
    figure, axes = plt.subplots(figsize=(12,10))
    basemap = Basemap(projection='robin', lon_0=0, lat_0=0, resolution='h', area_thresh=1000, ax=axes) 
    basemap.drawcoastlines(linewidth=0.5) 
    basemap.drawparallels(np.arange(-90.,120.,30.), labels=[1,0,0,0])
    basemap.drawmeridians(np.arange(0.,420.,60.), labels=[0,0,0,1])
    return figure, axes, basemap

def plot_mesh_data(interpolation, grid, basemap):
    colormesh = basemap.contourf(
        grid["x"], grid["y"], interpolation,
        levels=120, cmap='RdBu_r', antialiased=True
    )
    color_bar = basemap.colorbar(colormesh, location='bottom', pad="10%")
    return colormesh

def mask_greenland(axes, basemap, shp_path):
    basemap.readshapefile(shp_path, 'greenland', drawbounds=False)
    patches = []
    for shape in basemap.greenland:
        patches.append(Polygon(np.array(shape), True))
    mask = PatchCollection(patches, facecolor='white', edgecolor='k', linewidths=0.5, zorder=3)
    axes.add_collection(mask)

if __name__ == "__main__":
    df = load_data()
    base_data = get_data(df)
    # 加权南极观测点
    base_data = weight_antarctica_data(base_data, weight=4)
    figure, axes, basemap = prepare_map_plot()
    grid = generate_grid(base_data, basemap, delta=0.5)
    extended_data = extend_data(base_data)
    interpolation, interpolation_error = interpolate(extended_data, grid)
    plot_mesh_data(interpolation, grid, basemap)
    # 替换为格陵兰shp文件的实际路径
    mask_greenland(axes, basemap, r"C:/Users/XXX/Desktop/Greenland/GRL_adm0")
    plt.show()

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