Python读取NC格式数据绘制水汽通量等值线和风场
; 一、知识点
计算水汽通量，用到了metpy包，是一个地球科学计算包，内置了很多气象用到的计算函数
小知识点：
1.用湿度计算比湿
2.单位的使用
3.常量的使用，这里涉及了重力加速度g
二、代码拆分
导入包
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import my_class
from metpy.units import units
import metpy.constants
import metpy.calc
读取数据
ds = xr.open_dataset('download_201306.nc')
lat = ds.latitude
lon = ds.longitude
time = ds.time
u = ds['u']
v = ds['v']
tem = ds['t']- 273.15
rh = ds['r']
注意：这里读取的数据是全部的格点数据，但是我们画图用不了这么多，所以对数据做分割。只需要一部分
数据筛选分割
lat_range = lat[(lat>=22) & (lat33)]
lon_range = lon[(lon>=106) & (lon124)]
u_region = u.sel(longitude=lon_range, latitude=lat_range,time = '2013-06-28T00:00:00.000000000')
v_region = v.sel(longitude=lon_range, latitude=lat_range,time = '2013-06-28T00:00:00.000000000')
tem_region = tem.sel(longitude=lon_range, latitude=lat_range,time = '2013-06-28T00:00:00.000000000')
rh_region = rh.sel(longitude=lon_range, latitude=lat_range,time = '2013-06-28T00:00:00.000000000')
上面获得了世界时2013年6月28日00时的UV风场、温度、湿度
这里解释一下，作者下载的在分析资料里没有分层的数据， 只下载了850hPa一层数据，所以不需要设置高度。
数据计算
tem_region = tem_region * units.degC
pressure = 850 * units.hPa
dewpoint = metpy.calc.dewpoint_from_relative_humidity(tem_region, rh_region)
specific_humidity = metpy.calc.specific_humidity_from_dewpoint(pressure, dewpoint)*1000
qu = u_region * specific_humidity / metpy.constants.g
qv = v_region * specific_humidity / metpy.constants.g
a = np.sqrt(qu * qu + qv*qv)
这里的知识点：
1.温度和气压赋予单位，因为metpy的计算是有单位体系的
2.利用metpy.calc的内置函数计算物理量
dewpoint_from_relative_humidity：露点_from_相对湿度，参数是温度和湿度
specific_humidity_from_dewpoint：比湿_from_露点，参数是气压和露点。
3.利用水汽通量公式： 速度*比湿/重力加速度计算水汽通量的UV分量。
画图设置
fig = plt.figure(figsize=(9,6))
ax = fig.add_subplot(111,projection = ccrs.PlateCarree())
ax.set_extent([106, 124, 22, 33])
ax.set_xticks(np.arange(106,125,2),crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(22,34,1),crs=ccrs.PlateCarree())
xticks_str = ['106', '108', '110', '112', '114', '116', '118','120', '122', '124°E']
ax.set_xticklabels(xticks_str,fontsize = 11)
yticks_str = ['22 ','23 ','24 ','25 ','26 ','27 ','28 ','29 ','30 ','31 ','32 ','33°N']
ax.set_yticklabels(yticks_str,fontsize = 11)
my_class.readshapefile('bou2_4l.shp',linewidth=1,ax=ax)
画图设置这里就不多讲了，就是弄个画布、弄个子图，设置一下XY轴，加载个地图。
画图
ct=ax.contour(lon_range,lat_range,a,8,colors='k',linewidths=1,levels=range(0,28,2))
ax.clabel(ct,inline=True,fontsize=10,fmt='%.0f')
h1 = ax.quiver(lon_range[::4],lat_range[::4],qu[::4,::4],qv[::4,::4],
width = 0.002,
scale = 700,
pivot = 'tail'
)
ax.quiverkey(h1,
X=0.8, Y = -0.07,
U = 20,
angle = 0,
label='20',
labelpos='S',
color = 'k',labelcolor = 'k',
)
这里画了一个等值线，画了一个水汽通量的方向场，其实就是画风场的方法。
出图
plt.show()
三、完整代码
这里的完整代码是两张图，包含了一个循环，就是把两个时次的图放在一起。
'''
知识点：计算水汽通量，用到了metpy包，是一个地球科学计算包，内置了很多气象用到的计算函数
小知识点：1.用湿度计算比湿
2.单位的使用
3.常量的使用，这里涉及了重力加速度g
'''
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import my_class
from metpy.units import units
import metpy.constants
import metpy.calc
ds = xr.open_dataset('download_201306.nc')
lat = ds.latitude
lon = ds.longitude
time = ds.time
u = ds['u']
v = ds['v']
tem = ds['t']- 273.15
rh = ds['r']
lat_range = lat[(lat>=22) & (lat33)]
lon_range = lon[(lon>=106) & (lon124)]
fig = plt.figure(figsize=(18,9))
plt.subplots_adjust(left=0.07, right=0.90, top=0.95, bottom=0.05,wspace=0.2,hspace=0.1)
ax_a = fig.add_subplot(1,2,1,projection = ccrs.PlateCarree())
ax_b = fig.add_subplot(1,2,2,projection = ccrs.PlateCarree())
ax_list = [ax_a,ax_b]
ab = ['(a)','(b)']
i = 0
for time_i in time:
u_region = u.sel(longitude=lon_range, latitude=lat_range,time = time_i)
v_region = v.sel(longitude=lon_range, latitude=lat_range,time = time_i)
tem_region = tem.sel(longitude=lon_range, latitude=lat_range,time = time_i)
rh_region = rh.sel(longitude=lon_range, latitude=lat_range,time = time_i)
tem_region = tem_region * units.degC
pressure = 850 * units.hPa
dewpoint = metpy.calc.dewpoint_from_relative_humidity(tem_region, rh_region)
specific_humidity = metpy.calc.specific_humidity_from_dewpoint(pressure, dewpoint)*1000
qu = u_region * specific_humidity / metpy.constants.g
qv = v_region * specific_humidity / metpy.constants.g
a = np.sqrt(qu * qu + qv*qv)
ax_list[i].set_extent([106, 124, 22, 33])
ax_list[i].set_xticks(np.arange(106,125,2),crs=ccrs.PlateCarree())
ax_list[i].set_yticks(np.arange(22,34,1),crs=ccrs.PlateCarree())
xticks_str = ['106', '108', '110', '112', '114', '116', '118','120', '122', '124°E']
ax_list[i].set_xticklabels(xticks_str,fontsize = 11)
yticks_str = ['22 ','23 ','24 ','25 ','26 ','27 ','28 ','29 ','30 ','31 ','32 ','33°N']
ax_list[i].set_yticklabels(yticks_str,fontsize = 11)
my_class.readshapefile('bou2_4l.shp',linewidth=1,ax=ax_list[i])
ct=ax_list[i].contour(lon_range,lat_range,a,8,colors='k',linewidths=1,levels=range(0,28,2))
ax_list[i].clabel(ct,inline=True,fontsize=10,fmt='%.0f')
h1 = ax_list[i].quiver(lon_range[::4],lat_range[::4],qu[::4,::4],qv[::4,::4],
width = 0.002,
scale = 700,
pivot = 'tail'
)
ax_list[i].quiverkey(h1,
X=0.8, Y = -0.07,
U = 20,
angle = 0,
label='20',
labelpos='S',
color = 'k',labelcolor = 'k',
)
ax_list[i].text(0.05, 0.95, ab[i],transform=ax_list[i].transAxes, fontsize=11)
i += 1
plt.show()
Original: https://blog.csdn.net/weixin_42372313/article/details/125503942
Author: 野生的气象小流星
Title: Python读取NC格式数据绘制水汽通量等值线和和流场
原创文章受到原创版权保护。转载请注明出处：https://www.johngo689.com/730895/
转载文章受原作者版权保护。转载请注明原作者出处！