Python读取NC格式数据绘制水汽通量等值线和和流场

Python读取NC格式数据绘制水汽通量等值线和风场

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']

注:这里读取的数据是所有的栅格数据,但我们在绘图中不能用那么多,所以我们对数据进行了划分。只需要其中的一部分。

[En]

Note: the data read here is all the grid data, but we can’t use so much in drawing, so we divide the data. Only part of it is needed.

数据筛选分割


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',
             )

这里我们画了一条等值线,画出了水汽通量的方向场,这实际上是绘制风场的方法。

[En]

Here we draw an isoline and draw a direction field of water vapor flux, which is actually the way to draw the wind field.

出图

plt.show()

Python读取NC格式数据绘制水汽通量等值线和和流场

三、完整代码

这里的完整代码是两个图表,包含一个循环,即将这两个时间的图表放在一起。

[En]

The complete code here is two diagrams, containing a loop, that is, putting the diagrams of the two times together.

'''
知识点:计算水汽通量,用到了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/30674/

转载文章受原作者版权保护。转载请注明原作者出处!

(0)

大家都在看

发表回复

登录后才能评论
免费咨询
免费咨询
扫码关注
扫码关注
联系站长

站长Johngo!

大数据和算法重度研究者!

持续产出大数据、算法、LeetCode干货,以及业界好资源!

2022012703491714

微信来撩,免费咨询:xiaozhu_tec

分享本页
返回顶部