说到插值大概会提到日常用的scipy的linear和cubic,克里金插值等等 meteva也有插值功能,不论是站点插网格,网格插站点,还是网格插网格统统都有 本文主要测试meteva的IDW与cressman站点插值 并基于插值后的数据测试插值后的白化效果
版本:python3.9
%matplotlib inline
%load_ext autoreload
%autoreload 2
import meteva.base as meb
from cnmaps import get_adm_maps, draw_maps, clip_contours_by_map
import metpy.calc as mpcalc
import metpy.plots as mpplots
from metpy.units import units
from metpy.calc import reduce_point_density
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.io.shapereader import BasicReader
from cartopy.mpl.gridliner import LATITUDE_FORMATTER, LONGITUDE_FORMATTER
import cmaps
plt.rcParams['font.sans-serif'] = ['Source Han Sans CN']
plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号
读取数据并绘制站点图
# 读取Micaps站点数据
provinces = BasicReader('/home/mw/input/china1656/china_map/china_map/China_Province_2022.shp')
filename = "/home/mw/input/meteva2260/20210719120000.000" # 替换为你的micaps文件路径
sta = meb.read_stadata_from_micaps3(filename)
sta.head()
lons = sta["lon"].values
lats = sta["lat"].values
pr = sta["data0"].values * units.degC
# 创建图形
fig = plt.figure(figsize=(15, 12))
ax = fig.add_subplot(111,projection=ccrs.PlateCarree())
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.COASTLINE)
ax.set_extent([80, 130,20, 60])
# 创建站点图
stationplot = mpplots.StationPlot(ax, lons, lats, transform=ccrs.PlateCarree(),
fontsize=10)
# 绘制分布
stationplot.plot_parameter('NW',pr, formatter=lambda v: format(v, '.1f'))
ax.add_geometries(provinces.geometries(), linewidth=.5, edgecolor='black', crs=ccrs.PlateCarree(),
facecolor='none')
plt.show()
## 插值前要设置格点
grid1 = meb.grid([80,130,0.5],[20,60,0.5])
sta1 = meb.interp_sg_idw(sta,grid1,nearNum = 2)
sta2 = meb.interp_sg_cressman(sta,grid = grid1,r_list = [1000,200,100,50],nearNum = 100)
map_extend = [80, 130, 20, 60]
axs = meb.creat_axs(2, map_extend,ncol=2,sup_fontsize=7)
cmap = cmaps.radar_1
image = meb.add_contourf(axs[0], sta1,cmap =cmap,add_colorbar=False)
#image = meb.add_scatter_text(axs[0], sta, tag=0, font_size=5)
image = meb.add_contourf(axs[1], sta2,cmap =cmap)
#image = meb.add_scatter_text(axs[1], sta, tag=0, font_size=5)
左边为IDW,右边为cressman
map_extend = [110, 125, 35, 45]
axs = meb.creat_axs(2, map_extend,ncol=2,sup_fontsize=7)
image = meb.add_contourf(axs[0], sta1,cmap =cmap,add_colorbar=False,clip =["beijing","tianjin","hebei"])
#image = meb.add_scatter_text(axs[0], sta, tag=0, font_size=5)
image = meb.add_contourf(axs[1], sta2,cmap =cmap,clip =["beijing","tianjin","hebei"])
#image = meb.add_scatter_text(axs[1], sta, tag=0, font_size=5)
beijing = get_adm_maps(province='北京市', only_polygon=True, record='first')
tianjin = get_adm_maps(province='天津市', only_polygon=True, record='first')
hebei = get_adm_maps(province='河北省', only_polygon=True, record='first')
jingjinji = beijing + tianjin + hebei
# 创建图形
fig = plt.figure(figsize=(20, 12))
# 添加第一幅子图
ax1 = fig.add_subplot(221, projection=ccrs.PlateCarree())
ax1.set_extent([110, 125, 35, 45])
cmap1 = cmaps.radar
# 绘制填色分布图
im1 = ax1.contourf(sta1.lon,sta1.lat,sta1.values[0, 0,0,0,:,:],
cmap=cmap1,
transform=ccrs.PlateCarree(),
level=np.arange(0,28,2))
clip_contours_by_map(im1, jingjinji)
draw_maps(get_adm_maps(level='省'), linewidth=0.8, color='k')
fig.colorbar(im1, orientation='vertical')
# 添加第二幅子图
ax2 = fig.add_subplot(222, projection=ccrs.PlateCarree())
ax2.set_extent([110, 125, 35, 45])
# 绘制填色分布图
im2 = ax2.contourf(sta2.lon, sta2.lat, sta2.values[0, 0, 0, 0, :, :],
cmap=cmap1,
transform=ccrs.PlateCarree(),
level=np.arange(0, 28, 2))
# 添加色标
fig.colorbar(im2, orientation='vertical')
# 白化
clip_contours_by_map(im2, jingjinji)
draw_maps(get_adm_maps(level='省'), linewidth=0.8, color='k')
# 显示图形
plt.show()
meteva读取或插值后的数据为六维数组,因此绘图前需处理。metava好处在于便利,但需要注意其函数内部要素的排列比较严格,如报错时使用help(函数名)加以确认