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

用户11172986

发布于 2024-06-20 17:38:11

890

发布于 2024-06-20 17:38:11

文章被收录于专栏：气python风雨气python风雨

前言

gma库是洛大神写的一个地理库，他的公众号是洛的地理研学

小弟对这些开源大佬是非常感谢的，正是这种共享精神推动了python的发展.

言归正传，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]:

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

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

Out[47]:

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()