meteva站点插值填色与白化

用户11172986

发布于 2024-06-20 16:37:37

840

发布于 2024-06-20 16:37:37

前言

说到插值大概会提到日常用的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()

站点插值格点：IDW与cressman插值

IDW

## 插值前要设置格点
grid1 = meb.grid([80,130,0.5],[20,60,0.5])
sta1 = meb.interp_sg_idw(sta,grid1,nearNum = 2)

cressman插值

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

白化方法一：meteva内部函数白化（京津冀地区）

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)

白化方法二：基于cnmaps 白化（京津冀地区）

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