克里金法(Kriging)是依据协方差函数对随机过程/随机场进行空间建模和预测（插值）的回归算法。在特定的随机过程，例如固有平稳过程中，克里金法能够给出最优线性无偏估计（Best Linear Unbiased Prediction, BLUP），因此在地统计学中也被称为空间最优无偏估计器(spatial BLUP)。示意图如下：

像这种常用的算法一般都是有前人做过的，在github上一搜就能查到。对于克里金插值可以直接调用pykrige包进行Kriging插值计算。

代码

# coding: utf-8
import math
import numpy as np
import pandas as pd
import geopandas as gpd
import xarray as xr
import salem
import cartopy.crs as ccrs
import cartopy.feature as cfeat
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
from cartopy.io.shapereader import Reader, natural_earth
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
from pykrige.ok import OrdinaryKriging
import os

def create_map(fig):
    shp_path = r'/China_shp/'
    # --创建画图空间
    proj = ccrs.LambertConformal(central_longitude=110, central_latitude=35, 
                             standard_parallels=(30, 60))  # 创建坐标系

    ax = fig.subplots(1, 1, subplot_kw={'projection': proj})  
    # --设置shp
    provinces = cfeat.ShapelyFeature(
        Reader(shp_path + 'province.shp').geometries(),
        ccrs.PlateCarree(), 
        edgecolor='k',
        facecolor='none',
    )
    ax.add_feature(provinces, linewidth=0.6, zorder=2)
    # --设置范围
    ax.set_extent([80, 135, 15, 55], crs=ccrs.PlateCarree())
    # --关闭边框
    ax.spines['geo'].set_visible(False)
    # --设置刻度
    ax.tick_params(
        axis='both', 
        labelsize=5, 
        direction='out',
        labelbottom=False, 
        labeltop=False, 
        labelleft=False, 
        labelright=False
    )
    # --设置小地图
    left, bottom, width, height = 0.675, 0.13, 0.23, 0.27
    ax2 = fig.add_axes(
        [left, bottom, width, height], 
        projection=ccrs.PlateCarree()
    )
    provinces = cfeat.ShapelyFeature(
        Reader(shp_path + 'province.shp').geometries(),
        ccrs.PlateCarree(), 
        edgecolor='k',
        facecolor='gray',
    )
    ax2.add_feature(provinces, linewidth=0.6, zorder=2)
    ax2.set_extent([105, 125, 0, 25], crs=ccrs.PlateCarree())    
    return ax


def krige_interp(df):
    #插值成0.25*0.25
    grid_lon = np.arange(float(math.floor(df.lon.min())), float(math.ceil(df.lon.max())), 0.25)
    grid_lat = np.arange(float(math.floor(df.lat.min())), float(math.ceil(df.lat.max())), 0.25)
    lons = df["lon"].values
    lats = df["lat"].values
    data = df["data"].values
    
    #调用克里金插值函数
    OK = OrdinaryKriging(lons, lats, data,
                         variogram_model='gaussian',
                         coordinates_type="geographic",
                         nlags=6,)
    z, ss = OK.execute('grid', grid_lon, grid_lat)
    #转换成网格
    xgrid, ygrid = np.meshgrid(grid_lon, grid_lat)
    #将插值网格数据整理
    df_grid = pd.DataFrame(dict(lon=xgrid.flatten(), lat=ygrid.flatten()))
    #添加插值结果
    df_grid['data'] = z.flatten()
    df_grid = df_grid.set_index(['lat', 'lon'])
    ds = xr.Dataset.from_dataframe(df_grid).data
    return ds


def mask_region(ds, region='china'):
    shp_path = r'/China_shp'
    shp = gpd.read_file(shp_path + 'province.shp').set_index('省')
    if region != 'china':
        shp = shp.loc[region, :]
    ds_mask = ds.salem.roi(shape=shp)
    return ds_mask


def draw_fig(ds, levels, cmap, outpath='.'):
    fig = plt.figure(figsize=(10, 7), dpi=130)  # 创建页面
    
    im = ds.plot.contourf(
        ax=create_map(fig), 
        cmap=cmap, 
        transform=ccrs.PlateCarree(),
        levels=levels,
        add_colorbar=False
    )
    cbar_ax = fig.add_axes([0.13, 0.17, 0.3, 0.012])
    cb = fig.colorbar(im, cax=cbar_ax, orientation='horizontal')
    plt.savefig(f'{outpath}{os.sep}kriging_mask_chinamap.png')
    
if __name__ == "__main__":
    data_file = r'/testData.csv'
    levels = np.arange(116, 123, 0.5)
    cmap = 'Spectral_r'

    df = pd.read_csv(data_file, header=0, names=['lon', 'lat', 'data'])
    ds = krige_interp(df)
    ds_mask = mask_region(ds)
    draw_fig(ds_mask, levels, cmap,)

站点的密集程度对插值效果影响非常之大，由于中西部数据量比较少，所以插值结果也是没眼看。

也可以进一步调整，只留下部分省份，其余地区掩盖掉。

if __name__ == "__main__":
    data_file = r'./data.csv'
    levels = np.arange(116, 123, 0.5)
    cmap = 'Spectral_r'

    df = pd.read_csv(data_file, header=0, names=['lon', 'lat', 'data'])
    ds = krige_interp(df)
    ds_mask = mask_region(ds, region=['北京市', '天津市', '河北省', '山西省', '山东省', '辽宁省'])
    draw_fig(ds_mask, levels, cmap,)