两种降水站点数据克里金插值及可视化方法

前言

gma库是洛大神写的一个地理库，其中有许多可以使用的函数，今天简单介绍一下它克里金插值的使用，并与meteva进行对比

镜像：Python 3.9 GDAL3.4.3 核心函数：OrdinaryKriging

In [10]:

代码语言：javascript复制

pip install pykrige -i https://pypi.mirrors.ustc.edu.cn/simple/

代码语言：javascript复制

Looking in indexes: https://pypi.mirrors.ustc.edu.cn/simple/
Collecting pykrige
  Downloading https://mirrors.bfsu.edu.cn/pypi/web/packages/fa/5a/3bcb3ba5025e1047cb10867edb9659fcbbb2705fb57af2899a03a0f4195c/PyKrige-1.7.1-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (864 kB)
     |████████████████████████████████| 864 kB 2.2 MB/s eta 0:00:01
Requirement already satisfied: numpy<2,>=1.14.5 in /opt/conda/lib/python3.9/site-packages (from pykrige) (1.24.4)
Requirement already satisfied: scipy<2,>=1.1.0 in /opt/conda/lib/python3.9/site-packages (from pykrige) (1.10.1)
Installing collected packages: pykrige
Successfully installed pykrige-1.7.1
Note: you may need to restart the kernel to use updated packages.

数据读取

In [17]:

代码语言：javascript复制

%matplotlib inline
%load_ext autoreload
%autoreload 2
import meteva.base as meb
from pykrige.ok import OrdinaryKriging
import numpy as np
import xarray as xr

代码语言：javascript复制

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

In [2]:

代码语言：javascript复制

filename = "/home/mw/input/meteva2260/20210719120000.000"  # 替换为你的micaps文件路径
sta = meb.read_stadata_from_micaps3(filename)
sta.head()

代码语言：javascript复制

Out[2]:

	time	id	lon	lat	data0
0	2021-07-19 12:00:00	53763	112.1660	37.9044	0.3
1	2021-07-19 12:00:00	54279	126.5830	42.0500	0.1
2	2021-07-19 12:00:00	53770	112.3500	37.3567	0.1
3	2021-07-19 12:00:00	54287	128.0830	42.0167	0.5
4	2021-07-19 12:00:00	51730	81.2564	40.6064	1.3

METEVA

In [15]:

代码语言：javascript复制

## 插值前要设置格点
lons = sta['lon']
# 读取站点纬度
lats = sta['lat']
# 降水
data = sta['data0']
x = np.arange(105,125,0.1)
y= np.arange(25,50,0.1)

代码语言：javascript复制

In [26]:

代码语言：javascript复制

#构建测试数据 
OK = OrdinaryKriging(lons, lats, data, variogram_model='gaussian',nlags=6)
da0, ss = OK.execute('grid', x, y) 
da01 = xr.DataArray(da0, coords=[y,x], dims=['lat','lon']) 
print(da01)  #da0是一个DataArray格式数据
grd =meb.xarray_to_griddata(da01)

代码语言：javascript复制

<xarray.DataArray (lat: 250, lon: 200)>
array([[0.34398246, 0.22104452, 0.1802794 , ..., 2.00354938, 2.00354938,
        2.00354938],
       [0.34114398, 0.18142633, 0.11777748, ..., 2.00354938, 2.00354938,
        2.00354938],
       [0.40715308, 0.20817193, 0.11747919, ..., 2.00354938, 2.00354938,
        2.00354938],
       ...,
       [2.00354938, 2.00354938, 2.00354938, ..., 2.06695758, 2.03128827,
        2.01022228],
       [2.00354938, 2.00354938, 2.00354938, ..., 1.99143638, 1.9761684 ,
        1.97127591],
       [2.00354938, 2.00354938, 2.00354938, ..., 1.91555218, 1.92051165,
        1.93176057]])
Coordinates:
  * lat      (lat) float64 25.0 25.1 25.2 25.3 25.4 ... 49.5 49.6 49.7 49.8 49.9
  * lon      (lon) float64 105.0 105.1 105.2 105.3 ... 124.6 124.7 124.8 124.9

In [28]:

代码语言：javascript复制

map_extend = [105, 125, 25, 50]
axs = meb.creat_axs(1, map_extend,ncol=1,sup_fontsize=7)
image = meb.add_mesh(axs[0], grd ,add_colorbar=True)

代码语言：javascript复制

GMA

In [3]:

代码语言：javascript复制

import gma
from gma import io
from gma.smc import Interpolate
from gma.map import  plot, inres

Points = sta.loc[:, ['lon','lat']].values
Values = sta.loc[:, ['data0']].values

# 步骤1：反距离权重插值
KD = gma.smc.Interpolate.Kriging(Points, Values, Resolution = 0.1, 
                                 VariogramModel = 'Spherical', 
                                 VariogramParameters = None,
                                 KMethod = 'Ordinary',
                                 InProjection = 'EPSG:4326')
# 步骤2：将插值结果转换为 DataSet 数据集
KDDataSet = io.ReadArrayAsDataSet(KD.Data, Projection = 'WGS84', Transform = KD.Transform)

代码语言：javascript复制

In [5]:

代码语言：javascript复制

# 1.初始化一个地图框，并配置视图范围 
MapF = plot.MapFrame(Axes = None, Extent = [105, 25, 125, 50])

# 2.将内置的世界矢量图层添加到地图框
MapL1 = MapF.AddLayer(inres.WorldLayer.Country, FaceColor = 'none', LineWidth = 0.2, EdgeColor = 'black', Zorder = 2)
MapL2 = MapF.AddLayer(inres.WorldLayer.Ocean, FaceColor = '#BEE8FF', EdgeColor = 'none')

MapD1 = MapF.AddDataSetClassify(KDDataSet,
                                CMap = 'rainbow' )
# 3.获取经纬网
GridLines = MapF.AddGridLines(LONRange = (105, 130, 5), LATRange = (20, 60, 5))

# 4.添加地图整饰要素
AddCompass = MapF.AddCompass(LOC = (0.1, 0.8), Color = 'black')
ScaleBar = MapF.AddScaleBar(LOC = (0.1, 0.02), Width = 0.3, Color = 'black', FontSize = 7)

# 5.设置地图框边框
Frame = MapF.SetFrame()

代码语言：javascript复制

这里展现了如何创建xarray数组以及将xarray数组转为meteva可以可视化的griddata格式学习了这个即可实现快速可视化言归正传，两者大值分布仍然一致，但分辨率过高的gma低值分布明显不自然，相较上期IDW插值效果而言完结撒花

函数可视化配置数据数组

0 人点赞