使用python自动化绘制WRF模式输出的风场以及辐射

本脚本主要用来自动化处理WRF模式数据，可以根据自己指定的时间范围以及时间步长绘制相应的数据

1 导入库

import cmaps
import numpy as np
import glob
from netCDF4 import Dataset
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
import cartopy.crs as crs
from cartopy.feature import NaturalEarthFeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from wrf import (to_np, getvar, smooth2d, get_cartopy, cartopy_xlim,cartopy_ylim, latlon_coords,extract_global_attrs)
import proplot as pplt
import pandas as pd
import datetime
import os

2 选择数据范围

需要自己手动选择：初始时间、结束时间、时间步长

输入格式为：年-月-日-时、年-月-日-时、时
2022071000
2022071012
03

date = input()
frst = input()
step = input()path = r'G:/'+ date #只能是已经存在的文件目录且有数据才可以进行读取start = datetime.datetime.strptime(date,'%Y%m%d%H').strftime("%Y-%m-%d_%H_%M_%S")end = (datetime.datetime.strptime(date,'%Y%m%d%H')+datetime.timedelta(hours=int(frst))).strftime("%Y-%m-%d_%H_%M_%S")intp = (datetime.datetime.strptime(date,'%Y%m%d%H')+datetime.timedelta(hours=int(step))).strftime("%Y-%m-%d_%H_%M_%S")fstart = path+'/*'+start
fintp  = path+'/*'+intp
fend   = path+'/*'+endfile = path+'/*'filestart = glob.glob(fstart)
fileintp = glob.glob(fintp)
fileend = glob.glob(fend)filelist  = glob.glob(file)filelist.sort()   rstart = np.array(np.where(np.array(filelist)==filestart))[0][0]
rintp = np.array(np.where(np.array(filelist)==fileintp))[0][0]
rend   = np.array(np.where(np.array(filelist)==fileend))[0][0]

3定义读取数据函数以及绘图函数

定义函数读取风速资料U、V，以及辐射（辐射包括三部分）

将读取的资料进行绘图，并将绘图后的结果自动保存到指定路径下，方面后续绘制动图

def plot(i):ncfile = Dataset(i)#get u10 v10u10 = getvar(ncfile,'U10')v10 = getvar(ncfile,'V10')w10 = np.sqrt(u10*u10+v10*v10)#get net_radswdown = getvar(ncfile,'SWDOWN')glw    = getvar(ncfile,'GLW')olr    = getvar(ncfile,'OLR')rad    = swdown+glw+olrWE = extract_global_attrs(ncfile,'WEST-EAST_GRID_DIMENSION')['WEST-EAST_GRID_DIMENSION']SN = extract_global_attrs(ncfile,'SOUTH-NORTH_GRID_DIMENSION')['SOUTH-NORTH_GRID_DIMENSION']lev= extract_global_attrs(ncfile,'BOTTOM-TOP_GRID_DIMENSION')['BOTTOM-TOP_GRID_DIMENSION']t =u10.Time# time = pd.to_datetime((t)).strftime('%Y-%m-%d-%H')# Get the latitude and longitude pointslats, lons = latlon_coords(u10)# Create a figuref = pplt.figure(figsize=(5,5))ax = f.subplot(111, proj='lcyl')m=ax.contourf(to_np(lons), to_np(lats), to_np(w10),cmap=cmaps.MPL_RdYlBu_r,levels=np.arange(4,10.5,0.5),extend='both')ax.contour(to_np(lons), to_np(lats), to_np(rad),levels=8,linewidth=0.9,color='grey',# cmap=cmaps.MPL_RdYlBu_r,)step=4ax.quiver(to_np(lons)[::step,::step],to_np(lats)[::step,::step],to_np(u10)[::step,::step],to_np(v10)[::step,::step],)ax.format(ltitle='Sea surface WindSpeed (m/s)',coast=True,labels=True,lonlim=(120, 130),latlines=2,lonlines=2,latlim=(25, 35),)props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)bbox_props = dict(boxstyle="round", fc="w", ec="0.5", alpha=1)textstr = '\n'.join((r'$  Init:$'+ ncfile.START_DATE,r'$Valid:$' + i[-19:],# r'$\sigma=%.2f$' % (2, ))))ax.text(0.0, 1.2, textstr, transform=ax.transAxes, fontsize=9,verticalalignment='top', )ax.text(0.65, 1.2, 'U at 10 M(m s$^{-1}$)'+'\n'+'net_radiation(W m$^{-2}$)', transform=ax.transAxes, fontsize=9,verticalalignment='top', )ax.text(0.00, -0.07, ncfile.TITLE+'\n'+r'WE='+str(WE)+';SN='+str(SN)+';Lev='+str(lev)+';dx='+str(int(ncfile.DX))+'km;'+'Phys opt='+str(ncfile.MP_PHYSICS)+';PBL Opt='+str(ncfile.BL_PBL_PHYSICS)+';Cu Opt='+str(ncfile.CU_PHYSICS),fontsize=8,horizontalalignment='left',verticalalignment='top',transform=ax.transAxes,)f.colorbar(m, loc='r',label='',length=0.75,width=0.15)plt.show()savepath = r'J:/'+ dateif not os.path.exists(savepath):os.makedirs(savepath)f.savefig(savepath+'\\uv10'+i[-19:]+'.jpg')return

4 定义代码运行函数。

def plot_uv(fn):for i  in fn:plot(i)print('finish'+i)returnif __name__ == "__main__":plot_uv(fn)print('uvplot is finish')

5 个例测试

方面数据测试，取削封装函数运行。

path = r'G:\2022071000/*'file = glob.glob(path)file.sort()ncfile = Dataset(file[0])#get u10 v10
u10 = getvar(ncfile,'U10')
v10 = getvar(ncfile,'V10')
w10 = np.sqrt(u10*u10+v10*v10)
#get net_rad
swdown = getvar(ncfile,'SWDOWN')
glw    = getvar(ncfile,'GLW')
olr    = getvar(ncfile,'OLR')rad    = swdown+glw+olrWE = extract_global_attrs(ncfile,'WEST-EAST_GRID_DIMENSION')['WEST-EAST_GRID_DIMENSION']
SN = extract_global_attrs(ncfile,'SOUTH-NORTH_GRID_DIMENSION')['SOUTH-NORTH_GRID_DIMENSION']
lev= extract_global_attrs(ncfile,'BOTTOM-TOP_GRID_DIMENSION')['BOTTOM-TOP_GRID_DIMENSION']t =datetime.datetime.strptime(u10.Time.values,'%Y%m%d%H').strftime("%Y-%m-%d_%H_%M_%S")
# time = pd.to_datetime((t)).strftime('%Y-%m-%d-%H')# Get the latitude and longitude points
lats, lons = latlon_coords(u10)# Create a figure
f = pplt.figure(figsize=(5,5))
ax = f.subplot(111, proj='lcyl')m=ax.contourf(to_np(lons), to_np(lats), to_np(w10),cmap=cmaps.MPL_RdYlBu_r,)
ax.contour(to_np(lons), to_np(lats), to_np(rad),levels=8,linewidth=0.9,color='grey',# cmap=cmaps.MPL_RdYlBu_r,)
step=4
ax.quiver(to_np(lons)[::step,::step],to_np(lats)[::step,::step],to_np(u10)[::step,::step],to_np(v10)[::step,::step],)ax.format(ltitle='Sea surface WindSpeed (m/s)',coast=True,labels=True,lonlim=(120, 130),latlines=2,lonlines=2,latlim=(25, 35),)props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)bbox_props = dict(boxstyle="round", fc="w", ec="0.5", alpha=1)textstr = '\n'.join((r'$  Init:$'+ ncfile.START_DATE,r'$Valid:$' + file[],# r'$\sigma=%.2f$' % (2, ))))
ax.text(0.0, 1.2, textstr, transform=ax.transAxes, fontsize=9,verticalalignment='top', )ax.text(0.65, 1.2, 'U at 10 M(m s$^{-1}$)'+'\n'+'net_radiation(W m$^{-2}$)', transform=ax.transAxes, fontsize=9,verticalalignment='top', )ax.text(0.00, -0.07, ncfile.TITLE+'\n'+r'WE='+str(WE)+';SN='+str(SN)+';Lev='+str(lev)+';dx='+str(int(ncfile.DX))+'km;'+'Phys opt='+str(ncfile.MP_PHYSICS)+';PBL Opt='+str(ncfile.BL_PBL_PHYSICS)+';Cu Opt='+str(ncfile.CU_PHYSICS),fontsize=8,horizontalalignment='left',verticalalignment='top',transform=ax.transAxes,)f.colorbar(m, loc='r',label='',length=0.75,width=0.15)
plt.show()
# f.savefig('J:\\20220710\\'+i[35:]+'.png')# for i in file:#     plot_p(i)#     print('finish'+i)# print(i)


# -*- coding: utf-8 -*-
"""
Created on %(date)s@author: %(jixianpu)sEmail : 211311040008@hhu.edu.cnintroduction : keep learning althongh walk slowly
"""
"""
introduction :绘制近地面的风场和辐射
"""
import cmaps
import numpy as np
import glob
from netCDF4 import Dataset
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
import cartopy.crs as crs
from cartopy.feature import NaturalEarthFeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from wrf import (to_np, getvar, smooth2d, get_cartopy, cartopy_xlim,cartopy_ylim, latlon_coords,extract_global_attrs)
import proplot as pplt
import pandas as pd
import datetime
import osdate = input()
frst = input()
step = input()path = r'G:/'+ date #只能是已经存在的文件目录且有数据才可以进行读取start = datetime.datetime.strptime(date,'%Y%m%d%H').strftime("%Y-%m-%d_%H_%M_%S")end = (datetime.datetime.strptime(date,'%Y%m%d%H')+datetime.timedelta(hours=int(frst))).strftime("%Y-%m-%d_%H_%M_%S")intp = (datetime.datetime.strptime(date,'%Y%m%d%H')+datetime.timedelta(hours=int(step))).strftime("%Y-%m-%d_%H_%M_%S")fstart = path+'/*'+start
fintp  = path+'/*'+intp
fend   = path+'/*'+endfile = path+'/*'filestart = glob.glob(fstart)
fileintp = glob.glob(fintp)
fileend = glob.glob(fend)filelist  = glob.glob(file)filelist.sort()   rstart = np.array(np.where(np.array(filelist)==filestart))[0][0]
rintp = np.array(np.where(np.array(filelist)==fileintp))[0][0]
rend   = np.array(np.where(np.array(filelist)==fileend))[0][0]fn = filelist[rstart:rend:rintp]
# fn = filelist[rstart:rend:rintp]def plot(i):ncfile = Dataset(i)#get u10 v10u10 = getvar(ncfile,'U10')v10 = getvar(ncfile,'V10')w10 = np.sqrt(u10*u10+v10*v10)#get net_radswdown = getvar(ncfile,'SWDOWN')glw    = getvar(ncfile,'GLW')olr    = getvar(ncfile,'OLR')rad    = swdown+glw+olrWE = extract_global_attrs(ncfile,'WEST-EAST_GRID_DIMENSION')['WEST-EAST_GRID_DIMENSION']SN = extract_global_attrs(ncfile,'SOUTH-NORTH_GRID_DIMENSION')['SOUTH-NORTH_GRID_DIMENSION']lev= extract_global_attrs(ncfile,'BOTTOM-TOP_GRID_DIMENSION')['BOTTOM-TOP_GRID_DIMENSION']t =u10.Time# time = pd.to_datetime((t)).strftime('%Y-%m-%d-%H')# Get the latitude and longitude pointslats, lons = latlon_coords(u10)# Create a figuref = pplt.figure(figsize=(5,5))ax = f.subplot(111, proj='lcyl')m=ax.contourf(to_np(lons), to_np(lats), to_np(w10),cmap=cmaps.MPL_RdYlBu_r,levels=np.arange(4,10.5,0.5),extend='both')ax.contour(to_np(lons), to_np(lats), to_np(rad),levels=8,linewidth=0.9,color='grey',# cmap=cmaps.MPL_RdYlBu_r,)step=4ax.quiver(to_np(lons)[::step,::step],to_np(lats)[::step,::step],to_np(u10)[::step,::step],to_np(v10)[::step,::step],)ax.format(ltitle='Sea surface WindSpeed (m/s)',coast=True,labels=True,lonlim=(120, 130),latlines=2,lonlines=2,latlim=(25, 35),)props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)bbox_props = dict(boxstyle="round", fc="w", ec="0.5", alpha=1)textstr = '\n'.join((r'$  Init:$'+ ncfile.START_DATE,r'$Valid:$' + i[-19:],# r'$\sigma=%.2f$' % (2, ))))ax.text(0.0, 1.2, textstr, transform=ax.transAxes, fontsize=9,verticalalignment='top', )ax.text(0.65, 1.2, 'U at 10 M(m s$^{-1}$)'+'\n'+'net_radiation(W m$^{-2}$)', transform=ax.transAxes, fontsize=9,verticalalignment='top', )ax.text(0.00, -0.07, ncfile.TITLE+'\n'+r'WE='+str(WE)+';SN='+str(SN)+';Lev='+str(lev)+';dx='+str(int(ncfile.DX))+'km;'+'Phys opt='+str(ncfile.MP_PHYSICS)+';PBL Opt='+str(ncfile.BL_PBL_PHYSICS)+';Cu Opt='+str(ncfile.CU_PHYSICS),fontsize=8,horizontalalignment='left',verticalalignment='top',transform=ax.transAxes,)f.colorbar(m, loc='r',label='',length=0.75,width=0.15)plt.show()savepath = r'J:/'+ dateif not os.path.exists(savepath):os.makedirs(savepath)f.savefig(savepath+'\\uv10'+i[-19:]+'.jpg')returndef plot_uv(fn):for i  in fn:plot(i)print('finish'+i)returnif __name__ == "__main__":plot_uv(fn)print('uvplot is finish')

绘制结果如下图所示：

python--绘制WRF模式近地面风场以及辐射相关推荐

WRF后处理：python cartopy绘制土地利用/土地分类图//python绘制WRF下垫面类型（以北极为例）
下垫面类型对于WRF的地表过程十分重要,而在我们研究WRF的地表过程之前,需要对输入的土地利用类型进行一些绘制,以便后续的修改. 土地利用分类图的特点主要在于色标的设置,而在Github中,已有人根据 ...
python数据分析实战：近地面臭氧浓度预测
文章目录 1.背景 2.数据处理原始数据中包含哪些信息? 我们需要什么形式的数据? 面对异常数据该如何处理? 3.模型训练模型选择交叉验证调参优化结果代码 1.背景在上海隔离期间闲来无 ...
【气象数值模式及其数据处理】WRF模式与Python融合
当今从事气象及其周边相关领域的人员,常会涉及气象数值模式及其数据处理,无论是作为业务预报的手段.还是作为科研工具,气象数值模式与高效前后处理语言是一件非常重要的技能. WRF作为中尺度气象数值模式的佼 ...
WRF模式与Python融合技术在多领域中的应用及精美绘图
当今从事气象及其周边相关领域的人员,常会涉及气象数值模式及其数据处理,无论是作为业务预报的手段.还是作为科研工具,掌握气象数值模式与高效前后处理语言是一件非常重要的技能.WRF作为中尺度气象数值模式的 ...
Python地球科学领域应用：python处理遥感数据、站点数据、遥感水文数据、气候变化数据、WRF模式数据后处理、运行生态模型
点击查看原文>>>Python地球科学[赠CMIP6月/日数据.全球VIPPHEN物候数据.ERA5-LAND陆面再分析数据.遥感降水数据] >>>高精度气象模拟软 ...
【案例实践】WRF-Python融合技术：WRF 模式前后处理、自动化运行、数据处理、可视化绘图
[查看原文]Python在WRF模型自动化运行及前后处理中实践技术应用当今从事气象及其周边相关领域的人员,常会涉及气象数值模式及其数据处理,无论是作为业务预报的手段.还是作为科研工具,掌握气象数值模 ...
WRF模式应用：天气预报、模拟分析观测气温、降水、风场、水汽和湿度、土地利用变化、土壤及近地层能量水分通量、土壤、水体、植被等相关气象变量
查看原文>>>高精度气象模拟软件WRF(Weather Research Forecasting)技术及案例应用目录区域气候模式理论知识梳理 Linux操作系统WRF模式系统实际 ...
WRF模式可以做什么？天气预报、模拟分析观测气温、降水、风场、水汽和湿度、土地利用变化、土壤及近地层能量水分通量、土壤、水体、植被等相关气象变量
原文>>> 高精度气象模拟软件WRF(Weather Research Forecasting)技术及案例应用气候是多个领域(生态.水资源.风资源及碳中和等问题)的主要驱动因素,合 ...
Python绘制气象风场
title: Python绘制气象风场 date: 2021-08-05 21:01:52 tag: 气象风场数据来源为欧洲中心再分析数据ERA5 hourly data on pressure le ...

python--绘制WRF模式近地面风场以及辐射

使用python自动化绘制WRF模式输出的风场以及辐射

1 导入库

2 选择数据范围

3定义读取数据函数以及绘图函数

4 定义代码运行函数。

5 个例测试

python--绘制WRF模式近地面风场以及辐射相关推荐

最新文章

热门文章