简介压缩包内含一个基于 C/MFC 的 Kriging 空间插值等值线绘图工程面向 GIS、地质勘探等领域需要掌握空间插值和等值线图绘制的学习者与开发者。代码实现了数据预处理、半方差函数分析、Kriging 权重求解、插值计算以及等值线图形渲染等环节配套有可编译的工程文件也附带可执行程序便于直接观察效果。压缩包为 RAR 格式共 84 个文件主要包括 C 源文件.h/.cpp、编译中间产物.obj/.pch、工程配置.dsw/.dsp、程序图标与图片资源以及测试数据、结果说明文档整体约 1.77MB。已有 1214 人学习/下载。通过阅读源码可以理解普通克里金与泛克里金模型的差异、半方差函数参数的拟合方式以及如何将插值结果绘制为等值线工程内文件类型和目录结构清晰尤其是头文件与实现文件分离的模块组织适合作为相关课程设计、空间数据分析或二次开发的参考起点。1. Kriging画等值线图先把散点数据变成连续面再让等值线自己说话搞气象、地质、土壤环境的人大多遇到过这个尴尬手里攥着几十个站点的降雨量或污染物浓度落在图上只有一个个点领导要的是一张能看出“哪块儿超标、哪块儿安全”的等值线图。散点图撑不住反距离加权画出来又全是同心圆。Kriging克里金不是这种“算术平均值”的思路它基于空间自相关性先拟合半变异函数再按无偏最优估计去插值画出来的等值线图既有高低起伏又带误差估计。这份资源就是干这个的——用Python实现克里金插值到等值线出图的完整流程附带可以直接替换的示例数据。适合需要把离散观测值网格化、出图、写报告或做污染评价的从业者省掉从零调参的折腾。2. 克里金插值原理与选型为什么等值线图先要拟合半变异函数2.1 克里金不是距离反加权核心是半变异函数等值线图的本质是一个网格化的连续曲面Kriging做的是“从已知点到未知点的最优线性无偏估计”。它和反距离加权IDW最大的区别在于权重不是简单地取距离的倒数而是由半变异函数描述的空间相关性决定。半变异函数反映了距离越近、相似度越高这一规律以及这种相关性衰减到零的距离——变程range。如果数据在120公里内还有相关性那么超过这个距离的点就基本不参与插值权重会自动趋于零而IDW给所有点都分配权重距离大的点也会影响面结果就是平滑有余、细节失真。在资源的核心代码里第一步就是绘制经验半变异函数散点图再拟合一个理论模型。初看这一步觉得可有可无实际上它决定了后续等值线图的形态。如果半变异函数拟合得很糟糕网格上的插值结果就有系统性偏差画出来的等值线要么整体漂移要么局部出现不该有的“牛眼”。所以说Kriging出图能不能用50%的功夫在插值之前。2.2 球状、指数、高斯、线性四个常用模型的取舍资源里默认内置了四个半变异函数模型它们的表达式和适用场景差别不小。球状模型spherical在变程处相关性恰好降为零数学性质干净是最常用的通用模型适合土壤性质、降雨量这类相关距离明显的变量。指数模型exponential相关性渐近地趋近于零更平滑但在变程附近会有较好连续性适合地形高程、温度这类自相关衰减较慢的数据。高斯模型gaussian在原点附近特别平缓适合非常光滑的物理场但如果数据有噪声会过度拟合。线性模型linear没有固定变程简单但糙往往在数据量少、看不出明显变程时兜底用。选型没有绝对对错我的习惯是先用球状模型跑一遍再看交叉验证的误差如果残差有结构换指数或高斯对比。资源里提供了拟合报表直接比较AIC或RMS不用靠肉眼猜。2.3 网格分辨率、变程和块金值三个影响结果的参数很多人拿到代码后只改数据路径网格间距沿用默认结果出的图要么锯齿严重要么计算慢得离谱。网格分辨率要跟变程匹配一般取变程的1/101/20作为网格间距。如果变程是12公里网格间距设1公里比较合理设成0.1公里数据点之间全是外推网格数暴涨图也不见得更准。另一个参数是块金值nugget表示测量误差或微观变异。设为0时插值曲面会强行穿过每个观测点等值线图上出现“麻点”设置稍大一些曲面会适度平滑等值线更干净。资源里的拟合工具可以自动估计块金值但如果数据噪声大手调一下更可控。# 伪代码展示网格间距和变程的关系 range_estimate 12.0 # 从半变异函数拟合得到的变程单位公里 grid_spacing range_estimate / 15 # 经验值网格间距≈变程/15 print(f推荐网格间距{grid_spacing:.2f} km)这里range_estimate需要从拟合结果里读不能拍脑袋。grid_spacing设太小内存占用会指数上升设太大等值线会失去细节。代码注释里强调过这个比例实际跑数据时值得先打印变程看一眼。3. 从散点到等值线图用Python把克里金流程跑通3.1 数据准备坐标转换和缺失值检查这是整个流程里最枯燥但最要命的一步。我见过太多人拿着经纬度直接做克里金插值结果在高纬度地区网格变形严重等值线图被拉成扁椭圆。如果数据范围在几个公里量级可以用高斯-克吕格投影把经纬度转成平面坐标如果研究区跨度超过几百公里要考虑分带或使用UTM。资源的数据准备脚本里内置了pyproj的转换函数直接传入EPSG代码即可。另外观测点里混入非数值或空值会导致半变异函数拟合直接报错或者返回全NaN。检查一条都不能省。常见做法是先筛除缺失值和明显异常值比如降雨量为负再输出站点密度图确认覆盖范围免得后面网格化时出现大片无数据区域。import pandas as pd import pyproj df pd.read_csv(rainfall_data.csv) df_clean df.dropna(subset[lon, lat, value]) # 经纬度转投影坐标以EPSG:32650为例UTM 50N transformer pyproj.Transformer.from_crs(EPSG:4326, EPSG:32650, always_xyTrue) x, y transformer.transform(df_clean[lon].values, df_clean[lat].values) df_clean[x] x / 1000.0 # 转为公里方便变程单位可读 df_clean[y] y / 1000.0dropna(subset[...])同时检查经纬度和值三列任何一列缺失都会丢弃transformer.transform返回的是米除以1000变成公里这样后面克里金变程的单位就是公里好跟空间尺度对应。坐标转换这一步是全流程里最容易被跳过的但等值线图是否变形全看它。3.2 网格生成和克里金插值核心代码块网格生成推荐使用numpy的meshgrid范围取数据点的最小外接矩形再向外扩展少许缓冲。扩展量一般取变程的10%20%太多会导致外推区域忽悠人太少则图被裁剪到边缘。资源里有一个函数封装了网格生成并默认向外扩10%变程我觉得这个值在多数场景下很稳健。插值本身用pykrige的OrdinaryKriging。这里要传入观测点的坐标、值以及半变异函数模型。如果数据量大建议启用n_closest限制每个点只取邻近的20个观测点既加快计算又避免远处付作用。import numpy as np from pykrige.ok import OrdinaryKriging # 输入df_clean 已包含 x, y, value单位已转为公里 grid_x np.arange(west, east, grid_spacing) grid_y np.arange(south, north, grid_spacing) ok OrdinaryKriging( xdf_clean[x].values, ydf_clean[y].values, zdf_clean[value].values, variogram_modelspherical, nlags15, weightFalse, ) z_grid, ss_grid ok.execute(grid, grid_x, grid_y, n_closest20)nlags15表示经验半变异函数最多分15个距离段超出这个数拟合曲线会过于抖动n_closest20是克里金权重计算时只考虑最近的20个点避免了全数据量矩阵求逆计算速度提升明显。ss_grid是每个网格点的估计方差画不确定性等值线图时也可以拿来用这是克里金相对其他方法的一个额外好处。3.3 等值线绘制与渲染参数得到z_grid以后直接用matplotlib的contourf填充等值面再叠加contour画等值线。两个函数要分开配参数contourf控制色彩填充contour控制线宽和标注。最重要的参数是levels也就是等值线层级。很多人直接设levnp.linspace(z_min, z_max, 10)结果数据集中在低值段高值段只画出一两条超长的线图面很空。更好的做法是先看数据的分位数再决定levels。import matplotlib.pyplot as plt import matplotlib.ticker as ticker levels np.percentile(df_clean[value].values, np.linspace(10, 90, 9)) levels np.unique(np.concatenate(([df_clean[value].min()], levels, [df_clean[value].max()]))) fig, ax plt.subplots(figsize(8, 6)) cf ax.contourf(grid_x, grid_y, z_grid, levelslevels, cmapYlGnBu, alpha0.85) cs ax.contour(grid_x, grid_y, z_grid, levelslevels, colorsk, linewidths0.5) ax.clabel(cs, fmt%.1f, fontsize8) cbar fig.colorbar(cf, axax)np.percentile让层级按数据分布密度分布低值区和高值区都能被图层覆盖不会出现色块挤在一起的情况。clabel自动标注等值线数值如果图太密可以只标每第2条线避免标注叠在一起。最后加上站点散点图和边界就是一张能放进报告里的图。3.4 出图样式从能出图到出好图很多初学者的etc图能跑出来但受众是评审、业主或编辑美观度会影响可信度。资源和代码里给了几种样式模板一是把海岸线或行政边界叠加到等值线图上防止数据点落在境外或水域二是用基底地图basemap或cartopy显示背景但注意cartopy版本兼容问题三是当数据量级变化很大时把色标改成对数映射。from matplotlib.colors import LogNorm if data_range_ratio 20: norm LogNorm(vminz_min, vmaxz_max) else: norm None cf ax.contourf(grid_x, grid_y, z_grid, levelslevels, cmapYlGnBu, normnorm)data_range_ratio是最大值除以最小值。如果比值超过20线性色标会让低值区域完全看不出差异对数色标才合理。这个判断虽简单但能省掉很多“图出来没法看”的返工。4. 避坑指南克里金等值线图常见的五个翻车点4.1 半变异函数拟合失败数据不上正态变换现象半变异函数散点图乱得像云拟合曲线总是偏出来的等值线图一片色块之间没有过渡看起来像噪声云图。原因克里金本质假设数据近似正态分布偏态严重的原始数据比如降雨量右偏污染物浓度几个量级会让半变异函数的平方差被极大值主导相关性被掩埋。解决先做对数变换或Box-Cox变换插值完成后再反变换回原始量纲。资源里提供了自动判断偏度的函数偏度绝对值大于1就会提醒你先变换。4.2 经纬度未投影图被拉成“椭圆”现象站点纬度从北到南跨度10度经度跨度10度得到的等值线图格子不是方形而是横向或纵向拉长距离失真。原因直接把经纬度当平面坐标1度纬度和1度经度的物理距离在南北方不同克里金的“距离”概念就全乱了。解决参照3.1节用UTM投影把坐标转成米。注意跨分带时要用多带拼接工具否则带间隙位置会出现“缝合线”。4.3 网格范围随意设置边界内插外推失真现象等值线图在被边远站点包围的空白区域内出现大范围色块或者边界处等值线急剧扭曲。原因网格范围比数据范围大得多而克里金在数据覆盖范围之外没有观测点约束半变异函数会把远处变异全当成纯外推方差极大数值极为离谱。解决网格边界收缩到数据范围外扩变程的20%以内同时用ss_grid的估计方差给外推区域打上马赛克或标注不可信。等值线图只画在有数据约束的区域这是专业制图的基本素养。4.4 等值线层级不显示数据是离散型分布等值线被mask掉现象contourf画出的填充图正常但contour却只画出少量几条线或完全空白。原因如果z_grid里存在NaN常发生在数据覆盖区域外的网格点contour会默认忽略NaN区域并可能mask掉整条等值线另外如果levels设置不当比如层级数多于数据取值范围也会导致部分等值线不出现。解决先检查z_grid的NaN比例。如果NaN只出现在边缘用np.nanmax求覆盖范围并裁剪如果绞在中间说明搜索半径太小增大n_closest或变程上限。等值线层级数设置在515条比较稳妥过多会出现零碎短线。4.5 块金值设错热点都糊在一起现象等值线图在部分观测点周围出现一个个独立的圆形“小包”像得了麻疹或者反过来整个图过于平滑连已知的高值点都被抹平。原因块金值nugget设为零克里金严格插值穿过每个数据点局部噪声被当成真实特征放大块金值设置过大比如超过总方差的一半插值结果退化看向距离平均细节消失。解决观察半变异函数拟合截距。如果截距明显大于零但模型里写死nugget0就把nugget设为拟合出的截距值。资源里拟合函数会自动返回块金值不要忽略它直接改用这个值。5. 交叉验证与快速出图三个让克里金更靠谱的小技巧5.1 留一交叉验证选模型而不是拍脑袋工程上做克里金盲选模型的风险比想象中高。我会强制跑一遍留一交叉验证每次拿掉一个点用剩余的站点做插值再预测被拿掉的点累计所有误差。资源里提供了一条命令输出每个模型的MAE、RMSE和平均方差。我从那以后再也不敢说“球状模型一定好”了因为有一次高斯模型把RMSE降了20%。交叉验证的代码逻辑如下from pykrige.ok import OrdinaryKriging from sklearn.model_selection import LeaveOneOut import numpy as np def cv_score(x, y, z, model): rmses [] loo LeaveOneOut() for train_idx, test_idx in loo.split(x): ok OrdinaryKriging(x[train_idx], y[train_idx], z[train_idx], variogram_modelmodel) z_pred, _ ok.execute(points, x[test_idx], y[test_idx]) rmses.append((z[test_idx][0] - z_pred[0]) ** 2) return np.sqrt(np.mean(rmses)) for model in [spherical, exponential, gaussian]: print(model, cv_score(x, y, z, model))注意LeaveOneOut在数据量超过200个点时会很慢建议随机抽10%的点做交叉验证效果接近且速度翻倍。5.2 大数据量分块克里金别让内存爆掉当站点数超过5000个全数据矩阵的协方差求解会非常吃力普通笔记本直接卡死。我把研究区切成若干个有重叠的块每块单独插值最后用线性加权缝合。重叠区取变程的20%权重按到块中心的距离递减。这种方法不会损失太多精度但能让内存从爆掉变成平静运行。5.3 导出GeoTIFF方便GIS和Web端复用画完等值线图不只是交给甲方看很多场景下还需要放进ArcGIS或QGIS。这时可以把z_grid和网格坐标写进GeoTIFF保留地理投影信息这样同事直接拖进GIS就能出图。核心是给栅格重定义仿射变换参数保存为float32的栅格坐标。别小看这几个小技巧它们都是我从“等值线图交给别人后一问三不知”的状态里爬出来的教训。现在每完成一组克里金插值我都会留下交叉验证记录和半变异函数拟合报告这个文件既是自查依据也是给数据做“后悔药”。如果你也被散点数据搞到头大希望这几招能派上用场。希望帮到你。本文还有配套的精品资源点击获取
简介压缩包内含一个基于 C/MFC 的 Kriging 空间插值等值线绘图工程面向 GIS、地质勘探等领域需要掌握空间插值和等值线图绘制的学习者与开发者。代码实现了数据预处理、半方差函数分析、Kriging 权重求解、插值计算以及等值线图形渲染等环节配套有可编译的工程文件也附带可执行程序便于直接观察效果。压缩包为 RAR 格式共 84 个文件主要包括 C 源文件.h/.cpp、编译中间产物.obj/.pch、工程配置.dsw/.dsp、程序图标与图片资源以及测试数据、结果说明文档整体约 1.77MB。已有 1214 人学习/下载。通过阅读源码可以理解普通克里金与泛克里金模型的差异、半方差函数参数的拟合方式以及如何将插值结果绘制为等值线工程内文件类型和目录结构清晰尤其是头文件与实现文件分离的模块组织适合作为相关课程设计、空间数据分析或二次开发的参考起点。1. Kriging画等值线图先把散点数据变成连续面再让等值线自己说话搞气象、地质、土壤环境的人大多遇到过这个尴尬手里攥着几十个站点的降雨量或污染物浓度落在图上只有一个个点领导要的是一张能看出“哪块儿超标、哪块儿安全”的等值线图。散点图撑不住反距离加权画出来又全是同心圆。Kriging克里金不是这种“算术平均值”的思路它基于空间自相关性先拟合半变异函数再按无偏最优估计去插值画出来的等值线图既有高低起伏又带误差估计。这份资源就是干这个的——用Python实现克里金插值到等值线出图的完整流程附带可以直接替换的示例数据。适合需要把离散观测值网格化、出图、写报告或做污染评价的从业者省掉从零调参的折腾。2. 克里金插值原理与选型为什么等值线图先要拟合半变异函数2.1 克里金不是距离反加权核心是半变异函数等值线图的本质是一个网格化的连续曲面Kriging做的是“从已知点到未知点的最优线性无偏估计”。它和反距离加权IDW最大的区别在于权重不是简单地取距离的倒数而是由半变异函数描述的空间相关性决定。半变异函数反映了距离越近、相似度越高这一规律以及这种相关性衰减到零的距离——变程range。如果数据在120公里内还有相关性那么超过这个距离的点就基本不参与插值权重会自动趋于零而IDW给所有点都分配权重距离大的点也会影响面结果就是平滑有余、细节失真。在资源的核心代码里第一步就是绘制经验半变异函数散点图再拟合一个理论模型。初看这一步觉得可有可无实际上它决定了后续等值线图的形态。如果半变异函数拟合得很糟糕网格上的插值结果就有系统性偏差画出来的等值线要么整体漂移要么局部出现不该有的“牛眼”。所以说Kriging出图能不能用50%的功夫在插值之前。2.2 球状、指数、高斯、线性四个常用模型的取舍资源里默认内置了四个半变异函数模型它们的表达式和适用场景差别不小。球状模型spherical在变程处相关性恰好降为零数学性质干净是最常用的通用模型适合土壤性质、降雨量这类相关距离明显的变量。指数模型exponential相关性渐近地趋近于零更平滑但在变程附近会有较好连续性适合地形高程、温度这类自相关衰减较慢的数据。高斯模型gaussian在原点附近特别平缓适合非常光滑的物理场但如果数据有噪声会过度拟合。线性模型linear没有固定变程简单但糙往往在数据量少、看不出明显变程时兜底用。选型没有绝对对错我的习惯是先用球状模型跑一遍再看交叉验证的误差如果残差有结构换指数或高斯对比。资源里提供了拟合报表直接比较AIC或RMS不用靠肉眼猜。2.3 网格分辨率、变程和块金值三个影响结果的参数很多人拿到代码后只改数据路径网格间距沿用默认结果出的图要么锯齿严重要么计算慢得离谱。网格分辨率要跟变程匹配一般取变程的1/101/20作为网格间距。如果变程是12公里网格间距设1公里比较合理设成0.1公里数据点之间全是外推网格数暴涨图也不见得更准。另一个参数是块金值nugget表示测量误差或微观变异。设为0时插值曲面会强行穿过每个观测点等值线图上出现“麻点”设置稍大一些曲面会适度平滑等值线更干净。资源里的拟合工具可以自动估计块金值但如果数据噪声大手调一下更可控。# 伪代码展示网格间距和变程的关系 range_estimate 12.0 # 从半变异函数拟合得到的变程单位公里 grid_spacing range_estimate / 15 # 经验值网格间距≈变程/15 print(f推荐网格间距{grid_spacing:.2f} km)这里range_estimate需要从拟合结果里读不能拍脑袋。grid_spacing设太小内存占用会指数上升设太大等值线会失去细节。代码注释里强调过这个比例实际跑数据时值得先打印变程看一眼。3. 从散点到等值线图用Python把克里金流程跑通3.1 数据准备坐标转换和缺失值检查这是整个流程里最枯燥但最要命的一步。我见过太多人拿着经纬度直接做克里金插值结果在高纬度地区网格变形严重等值线图被拉成扁椭圆。如果数据范围在几个公里量级可以用高斯-克吕格投影把经纬度转成平面坐标如果研究区跨度超过几百公里要考虑分带或使用UTM。资源的数据准备脚本里内置了pyproj的转换函数直接传入EPSG代码即可。另外观测点里混入非数值或空值会导致半变异函数拟合直接报错或者返回全NaN。检查一条都不能省。常见做法是先筛除缺失值和明显异常值比如降雨量为负再输出站点密度图确认覆盖范围免得后面网格化时出现大片无数据区域。import pandas as pd import pyproj df pd.read_csv(rainfall_data.csv) df_clean df.dropna(subset[lon, lat, value]) # 经纬度转投影坐标以EPSG:32650为例UTM 50N transformer pyproj.Transformer.from_crs(EPSG:4326, EPSG:32650, always_xyTrue) x, y transformer.transform(df_clean[lon].values, df_clean[lat].values) df_clean[x] x / 1000.0 # 转为公里方便变程单位可读 df_clean[y] y / 1000.0dropna(subset[...])同时检查经纬度和值三列任何一列缺失都会丢弃transformer.transform返回的是米除以1000变成公里这样后面克里金变程的单位就是公里好跟空间尺度对应。坐标转换这一步是全流程里最容易被跳过的但等值线图是否变形全看它。3.2 网格生成和克里金插值核心代码块网格生成推荐使用numpy的meshgrid范围取数据点的最小外接矩形再向外扩展少许缓冲。扩展量一般取变程的10%20%太多会导致外推区域忽悠人太少则图被裁剪到边缘。资源里有一个函数封装了网格生成并默认向外扩10%变程我觉得这个值在多数场景下很稳健。插值本身用pykrige的OrdinaryKriging。这里要传入观测点的坐标、值以及半变异函数模型。如果数据量大建议启用n_closest限制每个点只取邻近的20个观测点既加快计算又避免远处付作用。import numpy as np from pykrige.ok import OrdinaryKriging # 输入df_clean 已包含 x, y, value单位已转为公里 grid_x np.arange(west, east, grid_spacing) grid_y np.arange(south, north, grid_spacing) ok OrdinaryKriging( xdf_clean[x].values, ydf_clean[y].values, zdf_clean[value].values, variogram_modelspherical, nlags15, weightFalse, ) z_grid, ss_grid ok.execute(grid, grid_x, grid_y, n_closest20)nlags15表示经验半变异函数最多分15个距离段超出这个数拟合曲线会过于抖动n_closest20是克里金权重计算时只考虑最近的20个点避免了全数据量矩阵求逆计算速度提升明显。ss_grid是每个网格点的估计方差画不确定性等值线图时也可以拿来用这是克里金相对其他方法的一个额外好处。3.3 等值线绘制与渲染参数得到z_grid以后直接用matplotlib的contourf填充等值面再叠加contour画等值线。两个函数要分开配参数contourf控制色彩填充contour控制线宽和标注。最重要的参数是levels也就是等值线层级。很多人直接设levnp.linspace(z_min, z_max, 10)结果数据集中在低值段高值段只画出一两条超长的线图面很空。更好的做法是先看数据的分位数再决定levels。import matplotlib.pyplot as plt import matplotlib.ticker as ticker levels np.percentile(df_clean[value].values, np.linspace(10, 90, 9)) levels np.unique(np.concatenate(([df_clean[value].min()], levels, [df_clean[value].max()]))) fig, ax plt.subplots(figsize(8, 6)) cf ax.contourf(grid_x, grid_y, z_grid, levelslevels, cmapYlGnBu, alpha0.85) cs ax.contour(grid_x, grid_y, z_grid, levelslevels, colorsk, linewidths0.5) ax.clabel(cs, fmt%.1f, fontsize8) cbar fig.colorbar(cf, axax)np.percentile让层级按数据分布密度分布低值区和高值区都能被图层覆盖不会出现色块挤在一起的情况。clabel自动标注等值线数值如果图太密可以只标每第2条线避免标注叠在一起。最后加上站点散点图和边界就是一张能放进报告里的图。3.4 出图样式从能出图到出好图很多初学者的etc图能跑出来但受众是评审、业主或编辑美观度会影响可信度。资源和代码里给了几种样式模板一是把海岸线或行政边界叠加到等值线图上防止数据点落在境外或水域二是用基底地图basemap或cartopy显示背景但注意cartopy版本兼容问题三是当数据量级变化很大时把色标改成对数映射。from matplotlib.colors import LogNorm if data_range_ratio 20: norm LogNorm(vminz_min, vmaxz_max) else: norm None cf ax.contourf(grid_x, grid_y, z_grid, levelslevels, cmapYlGnBu, normnorm)data_range_ratio是最大值除以最小值。如果比值超过20线性色标会让低值区域完全看不出差异对数色标才合理。这个判断虽简单但能省掉很多“图出来没法看”的返工。4. 避坑指南克里金等值线图常见的五个翻车点4.1 半变异函数拟合失败数据不上正态变换现象半变异函数散点图乱得像云拟合曲线总是偏出来的等值线图一片色块之间没有过渡看起来像噪声云图。原因克里金本质假设数据近似正态分布偏态严重的原始数据比如降雨量右偏污染物浓度几个量级会让半变异函数的平方差被极大值主导相关性被掩埋。解决先做对数变换或Box-Cox变换插值完成后再反变换回原始量纲。资源里提供了自动判断偏度的函数偏度绝对值大于1就会提醒你先变换。4.2 经纬度未投影图被拉成“椭圆”现象站点纬度从北到南跨度10度经度跨度10度得到的等值线图格子不是方形而是横向或纵向拉长距离失真。原因直接把经纬度当平面坐标1度纬度和1度经度的物理距离在南北方不同克里金的“距离”概念就全乱了。解决参照3.1节用UTM投影把坐标转成米。注意跨分带时要用多带拼接工具否则带间隙位置会出现“缝合线”。4.3 网格范围随意设置边界内插外推失真现象等值线图在被边远站点包围的空白区域内出现大范围色块或者边界处等值线急剧扭曲。原因网格范围比数据范围大得多而克里金在数据覆盖范围之外没有观测点约束半变异函数会把远处变异全当成纯外推方差极大数值极为离谱。解决网格边界收缩到数据范围外扩变程的20%以内同时用ss_grid的估计方差给外推区域打上马赛克或标注不可信。等值线图只画在有数据约束的区域这是专业制图的基本素养。4.4 等值线层级不显示数据是离散型分布等值线被mask掉现象contourf画出的填充图正常但contour却只画出少量几条线或完全空白。原因如果z_grid里存在NaN常发生在数据覆盖区域外的网格点contour会默认忽略NaN区域并可能mask掉整条等值线另外如果levels设置不当比如层级数多于数据取值范围也会导致部分等值线不出现。解决先检查z_grid的NaN比例。如果NaN只出现在边缘用np.nanmax求覆盖范围并裁剪如果绞在中间说明搜索半径太小增大n_closest或变程上限。等值线层级数设置在515条比较稳妥过多会出现零碎短线。4.5 块金值设错热点都糊在一起现象等值线图在部分观测点周围出现一个个独立的圆形“小包”像得了麻疹或者反过来整个图过于平滑连已知的高值点都被抹平。原因块金值nugget设为零克里金严格插值穿过每个数据点局部噪声被当成真实特征放大块金值设置过大比如超过总方差的一半插值结果退化看向距离平均细节消失。解决观察半变异函数拟合截距。如果截距明显大于零但模型里写死nugget0就把nugget设为拟合出的截距值。资源里拟合函数会自动返回块金值不要忽略它直接改用这个值。5. 交叉验证与快速出图三个让克里金更靠谱的小技巧5.1 留一交叉验证选模型而不是拍脑袋工程上做克里金盲选模型的风险比想象中高。我会强制跑一遍留一交叉验证每次拿掉一个点用剩余的站点做插值再预测被拿掉的点累计所有误差。资源里提供了一条命令输出每个模型的MAE、RMSE和平均方差。我从那以后再也不敢说“球状模型一定好”了因为有一次高斯模型把RMSE降了20%。交叉验证的代码逻辑如下from pykrige.ok import OrdinaryKriging from sklearn.model_selection import LeaveOneOut import numpy as np def cv_score(x, y, z, model): rmses [] loo LeaveOneOut() for train_idx, test_idx in loo.split(x): ok OrdinaryKriging(x[train_idx], y[train_idx], z[train_idx], variogram_modelmodel) z_pred, _ ok.execute(points, x[test_idx], y[test_idx]) rmses.append((z[test_idx][0] - z_pred[0]) ** 2) return np.sqrt(np.mean(rmses)) for model in [spherical, exponential, gaussian]: print(model, cv_score(x, y, z, model))注意LeaveOneOut在数据量超过200个点时会很慢建议随机抽10%的点做交叉验证效果接近且速度翻倍。5.2 大数据量分块克里金别让内存爆掉当站点数超过5000个全数据矩阵的协方差求解会非常吃力普通笔记本直接卡死。我把研究区切成若干个有重叠的块每块单独插值最后用线性加权缝合。重叠区取变程的20%权重按到块中心的距离递减。这种方法不会损失太多精度但能让内存从爆掉变成平静运行。5.3 导出GeoTIFF方便GIS和Web端复用画完等值线图不只是交给甲方看很多场景下还需要放进ArcGIS或QGIS。这时可以把z_grid和网格坐标写进GeoTIFF保留地理投影信息这样同事直接拖进GIS就能出图。核心是给栅格重定义仿射变换参数保存为float32的栅格坐标。别小看这几个小技巧它们都是我从“等值线图交给别人后一问三不知”的状态里爬出来的教训。现在每完成一组克里金插值我都会留下交叉验证记录和半变异函数拟合报告这个文件既是自查依据也是给数据做“后悔药”。如果你也被散点数据搞到头大希望这几招能派上用场。希望帮到你。本文还有配套的精品资源点击获取
