- 积分
- 1820
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2018-5-15
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
本帖最后由 Jesse123 于 2022-8-1 22:55 编辑
用wrfchem输出的物理量画图,兰伯特投影,数据和地图不对应(尤其是台湾),为啥呢???
代码如下:
mapdata=shpreader.Reader(r'E:\work\map\chinamap\bou2_4m\bou2_4p.shp')
f= xr.open_dataset(r'D:\A_meteo_workspace\wrfchem\geo_em.d02.nc') #ini(tional)为原始地形植被文件
lon = f['XLONG_M'][0,0,:]
lat = f['XLAT_M'][0,:,0]
file = xr.open_dataset(r'D:\A_meteo_workspace\wrfchem\wrfchemi_d02_2020-09-09_00-00-00.nc')
pm25i = np.squeeze(file['E_PM25I'][0,0,:,:])
pm25j = np.squeeze(file['E_PM25J'][0,0,:,:])
pm25 = pm25i+pm25j
#截取色标
level=[0,0.005,0.01,0.03,0.05,0.1,0.2,0.35,0.5,2.0]
cmap=cmaps.BlAqGrYeOrRe #获取Spectral色条,Spectral_r即为反色
newcolors=cm.get_cmap(cmap,10)
norm = mpl.colors.BoundaryNorm(level, cmap.N)
#开始绘图
plt.rc('font',family= 'Times New Roman')
# proj = ccrs.PlateCarree()
proj = ccrs.LambertConformal(central_longitude=114.5)
fig, ax1 = plt.subplots(figsize=(6.6, 6.6), subplot_kw=dict(projection=proj),dpi=400)
# ax1.set_xticks(np.arange(108.0, 121.57, 2),crs=ccrs.PlateCarree())
# ax1.set_yticks(np.arange(18.0, 28.3, 1),crs=ccrs.PlateCarree())
ax1.tick_params(axis='both',labelsize=12,direction='out',length=2.5,width=0.6,right=False,top=False)#修改刻度样式
ax1.set_extent([108.45,120.8,18.01,28.01],crs=ccrs.PlateCarree())
ax1.text(0.93,0.02,'d02',fontsize=16,transform=ax1.transAxes)
ax1.add_geometries(mapdata.geometries(),crs=ccrs.PlateCarree(),lw=0.9,facecolor='none',edgecolor='k',zorder=5)
a=ax1.contourf(lon, lat,pm25,levels=level,cmap=newcolors,norm=norm,extend='both',transform=ccrs.PlateCarree())
# 坐标与经纬网格(兰伯特投影)
fig.canvas.draw()
xticks = list(range(-180, 180, 2))
yticks = list(range(-90, 90, 2))
ax1.gridlines(xlocs=xticks, ylocs=yticks, linewidth=1.2, linestyle='--')
ax1.xaxis.set_major_formatter(LONGITUDE_FORMATTER)
ax1.yaxis.set_major_formatter(LATITUDE_FORMATTER)
lambert_ticks.lambert_xticks(ax1, xticks)
lambert_ticks.lambert_yticks(ax1, yticks)
plt.tick_params(labelsize=13)
#色标
a_level=[0,0.005,0.01,0.03,0.05,0.1,0.2,0.35,0.5,2.0]
cbposition=fig.add_axes([0.14, 0.12, 0.75, 0.013])
cb=plt.colorbar(a,ticks = a_level,orientation='horizontal',cax=cbposition,shrink = 0.8)
cb.set_label(r'$\mathregular{PM_{2.5}}/$'+r' '+r'$\mathregular{[ug/{m^3}*m/s]}$',fontdict={'weight':'normal','size':10})
cb.ax.tick_params(labelsize=10)
|
-
|