前言
gma库是洛大神写的一个地理库，
其中有许多可以使用的函数，
今天简单介绍一下它克里金插值的使用，并与meteva进行对比
镜像：Python 3.9 GDAL3.4.3
核心函数：OrdinaryKriging
In [10]:
pip install pykrige -i https://pypi.mirrors.ustc.edu.cn/simple/
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]:
%matplotlib inline%load_ext autoreload%autoreload 2import meteva.base as mebfrom pykrige.ok import OrdinaryKrigingimport numpy as npimport xarray as xr
The autoreload extension is already loaded. To reload it, use:
%reload_ext autoreload
In [2]:
filename = "/home/mw/input/meteva2260/20210719120000.000" # 替换为你的micaps文件路径sta = meb.read_stadata_from_micaps3(filename)sta.head()
Out[2]:
| level | time | dtime | id | lon | lat | data0 | |
|---|---|---|---|---|---|---|---|
| 0 | 0 | 2021-07-19 12:00:00 | 0 | 53763 | 112.1660 | 37.9044 | 0.3 |
| 1 | 0 | 2021-07-19 12:00:00 | 0 | 54279 | 126.5830 | 42.0500 | 0.1 |
| 2 | 0 | 2021-07-19 12:00:00 | 0 | 53770 | 112.3500 | 37.3567 | 0.1 |
| 3 | 0 | 2021-07-19 12:00:00 | 0 | 54287 | 128.0830 | 42.0167 | 0.5 |
| 4 | 0 | 2021-07-19 12:00:00 | 0 | 51730 | 81.2564 | 40.6064 | 1.3 |
METEVA
In [15]:
## 插值前要设置格点lons = sta['lon']# 读取站点纬度lats = sta['lat']# 降水data = sta['data0']x = np.arange(105,125,0.1)y= np.arange(25,50,0.1)
In [26]:
#构建测试数据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)
<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]:
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)
GMA
In [3]:
import gmafrom gma import iofrom gma.smc import Interpolatefrom gma.map import plot, inresPoints = sta.loc[:, ['lon','lat']].valuesValues = 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)
In [5]:
# 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()
这里展现了如何创建xarray数组以及将xarray数组转为meteva可以可视化的griddata格式
学习了这个即可实现快速可视化
言归正传，两者大值分布仍然一致，但分辨率过高的gma低值分布明显不自然，相较上期IDW插值效果而言
完结撒花
