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

前言

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

镜像：Python 3.9 GDAL3.4.3

数据读取

In [2]:

%matplotlib inline
%load_ext autoreload
%autoreload 2
import meteva.base as meb

In [3]:

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

Out[3]:

	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 [47]:

## 插值前要设置格点
grid1 = meb.grid([105,135,0.1],[25,55,0.1])
sta1 = meb.interp_sg_idw(sta,grid1,nearNum = 5)
sta1

Out[47]:

xarray.DataArray

'data0'

• member: 1
• level: 1
• time: 1
• dtime: 1
• lat: 301
• lon: 301

• array([[[[[[0.40511298, 0.39965576, 0.38122562, ..., 0.        ,
              0.        , 0.        ],
             [0.40206352, 0.39143518, 0.36779505, ..., 0.        ,
              0.        , 0.        ],
             [0.3975288 , 0.36514568, 0.33766782, ..., 0.        ,
              0.        , 0.        ],
             ...,
             [0.        , 0.        , 0.        , ..., 2.95151114,
              2.95273232, 2.95393586],
             [0.        , 0.        , 0.        , ..., 2.95258093,
              2.95379257, 2.95498729],
             [0.        , 0.        , 0.        , ..., 2.95366621,
              2.9548676 , 2.95605254]]]]]])

• Coordinates:

  • member

    (member)

    <U5

    'data0'

  • level

    (level)

    int64

    0

  • time

    (time)

    datetime64[ns]

    2021-07-19T12:00:00

  • dtime

    (dtime)

    int64

    0

  • lat

    (lat)

    float64

    25.0 25.1 25.2 ... 54.8 54.9 55.0

  • lon

    (lon)

    float64

    105.0 105.1 105.2 ... 134.9 135.0

• Indexes: (6)
• Attributes:

  dtime_units :

  hour

  data_start_columns :

  6

In [48]:

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

GMA

In [51]:

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：反距离权重插值
IDWD = Interpolate.IDW(Points, Values, Resolution = 0.1, InProjection = 'WGS84')
# 步骤2：将插值结果转换为 DataSet 数据集
IDWDataSet = io.ReadArrayAsDataSet(IDWD.Data, Projection = 'WGS84', Transform = IDWD.Transform)

In [52]:

# 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(IDWDataSet,
                                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()

当然就使用插值函数上两者都相对便利，而在可视化方面则是meteva对新手更优化，代码量二三行即可

gma在可视化方面则是需要设置更多的参数，估计需要一段时间上手之后才能画出更加美观的图

就插值结果而言两者分布是较为一致的，但是gma的resolution设为0.1和meteva设为0.1明显分辨率不一。还望有懂的前辈解惑一下。

还有则是个人的建议，希望gma插值后的数据格式可以选择为xarray，这样更加便利。