读取rmaps数据，并进行了等距离转等经纬度处理

lixiliang2008

f = addfile("E:/rmaps/rmapsnow_anal_init-20220727070000_3d.nc.gz")
x = f['X'][:]
y = f['Y'][:]
lat = 39.34 + y/30.0*0.27778
lon = 116.5 + x/(30.0*cos(40.95890833/180.0*3.1415926))*0.27778
lon, lat = meshgrid(lon, lat)
figure(figsize=[900, 630], newfig=False)
axesm(position=[0,0.53,0.54,0.46],tickfontsize=12,tickbold=True)

udata = f['uwind'][0,18,:,:]
vdata = f['vwind'][0,18,:,:]
speed = sqrt(udata * udata + vdata * vdata)
rdata = f['rh'][0,4,:,:]
levs = arange(60,100,10)
cols = [(255,255,255),(0,229,238),(0,238,0),(0,205,0),(0,160,0),(0,130,0)]
layer1 = barbsm(lon[::12,::12],lat[::12,::12],
         udata[::12,::12],vdata[::12,::12],speed[::12,::12],
         color='k',size = 9)
#plot map
mlayer = shaperead('F:\Meteinfo\MeteoInfo\map\cn_province.shp')
geoshow(mlayer, edgecolor='black', size=1.2)
xticks([],fontsize=12,bold=True)
yticks([36,37,38,39,40,41,42,],fontsize=12,bold=True)
ylim(lat.min(),lat.max())
xlim(lon.min(),lon.max())

f = addfile("E:/rmaps/rmapsnow_anal_init-20220727080000_3d.nc.gz")
x = f['X'][:]
y = f['Y'][:]
lat = 39.34 + y/30.0*0.27778
lon = 116.5 + x/(30.0*cos(40.95890833/180.0*3.1415926))*0.27778
lon, lat = meshgrid(lon, lat)

axesm(position=[0.4,0.53,0.54,0.46],tickfontsize=12,tickbold=True)

udata = f['uwind'][0,18,:,:]
vdata = f['vwind'][0,18,:,:]
speed = sqrt(udata * udata + vdata * vdata)
rdata = f['rh'][0,4,:,:]
levs = arange(60,100,10)
cols = [(255,255,255),(0,229,238),(0,238,0),(0,205,0),(0,160,0),(0,130,0)]
layer1 = barbsm(lon[::12,::12],lat[::12,::12],
         udata[::12,::12],vdata[::12,::12],speed[::12,::12],
         color='k',size = 9)
#plot map
mlayer = shaperead('F:\Meteinfo\MeteoInfo\map\cn_province.shp')
geoshow(mlayer, edgecolor='black', size=1.2)
xticks([],fontsize=12,bold=True)
yticks([],fontsize=12,bold=True)
ylim(lat.min(),lat.max())
xlim(lon.min(),lon.max())

f = addfile("E:/rmaps/rmapsnow_anal_init-20220727090000_3d.nc.gz")
x = f['X'][:]
y = f['Y'][:]
lat = 39.34 + y/30.0*0.27778
lon = 116.5 + x/(30.0*cos(40.95890833/180.0*3.1415926))*0.27778
lon, lat = meshgrid(lon, lat)

axesm(position=[0,0.04,0.54,0.46],tickfontsize=12,tickbold=True)

udata = f['uwind'][0,18,:,:]
vdata = f['vwind'][0,18,:,:]
speed = sqrt(udata * udata + vdata * vdata)
rdata = f['rh'][0,4,:,:]
levs = arange(60,100,10)
cols = [(255,255,255),(0,229,238),(0,238,0),(0,205,0),(0,160,0),(0,130,0)]
layer1 = barbsm(lon[::12,::12],lat[::12,::12],
         udata[::12,::12],vdata[::12,::12],speed[::12,::12],
         color='k',size = 9)
#plot map
mlayer = shaperead('F:\Meteinfo\MeteoInfo\map\cn_province.shp')
geoshow(mlayer, edgecolor='black', size=1.2)
xticks([113,114,115,116,117,118,119,120],fontsize=12,bold=True)
yticks([36,37,38,39,40,41,42,],fontsize=12,bold=True)
ylim(lat.min(),lat.max())
xlim(lon.min(),lon.max())

f = addfile("E:/rmaps/rmapsnow_anal_init-20220727100000_3d.nc.gz")
x = f['X'][:]
y = f['Y'][:]
lat = 39.34 + y/30.0*0.27778
lon = 116.5 + x/(30.0*cos(40.95890833/180.0*3.1415926))*0.27778
lon, lat = meshgrid(lon, lat)

axesm(position=[0.4,0.04,0.54,0.46],tickfontsize=12,tickbold=True)

udata = f['uwind'][0,18,:,:]
vdata = f['vwind'][0,18,:,:]
speed = sqrt(udata * udata + vdata * vdata)
rdata = f['rh'][0,4,:,:]
levs = arange(60,100,10)
cols = [(255,255,255),(0,229,238),(0,238,0),(0,205,0),(0,160,0),(0,130,0)]
layer1 = barbsm(lon[::12,::12],lat[::12,::12],
         udata[::12,::12],vdata[::12,::12],speed[::12,::12],
         color='k',size = 9)
#plot map
mlayer = shaperead('F:\Meteinfo\MeteoInfo\map\cn_province.shp')
geoshow(mlayer, edgecolor='black', size=1.2)
xticks([113,114,115,116,117,118,119,120],fontsize=12,bold=True)
yticks([],fontsize=12,bold=True)
ylim(lat.min(),lat.max())
xlim(lon.min(),lon.max())

MeteoInfo

rmaps数据是什么投影？addfile生成的 f 对象，f.proj 可以看看投影信息，如果能正确读取可以用MeteoInfoLab中的投影函数来进行坐标转换。此外绘图时对于MeteoInfo自带的shape文件不需要shaperead，只要在geoshow中给出shape文件名即可，比如 geoshow('cn_province', edgecolor='k', size=1.2)

lixiliang2008

应该是等经纬度数据，不带投影信息。