气象绘图—模仿绘制风速水汽通量图

沐雨淋觞

前言

本文是对气象业务网中，水汽通量图的模仿绘制。首发CSDN，后续气象相关绘图会持续更新。

一、读取数据

=数据来源=

ECMWF再分析数据

https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-pressure-levels?tab=form

目前CDS正在更新，暂无法下载（自动回复：请不要使用迅雷等下载工具，点我查看下载帮助）（自动回复：请不要使用迅雷等下载工具，点我查看下载帮助）。

“由于将数据存储的所有层（云、数据、软件）迁移到新的、最先进的硬件和软件基础设施上非常复杂，同时还要保持所有系统运行，每天提供超过 150
TB 的数据并支持 30 万用户，因此不可避免地会出现一些延迟，感谢您的耐心等待。”

1.1 读取netCDF文件

1. import xarray as xr
   
2. import numpy as np
   
3. filename = r'G:\document\document\biancheng\data\haikui_UVqVow.nc'
   
4. data     = xr.open_dataset(filename)

复制代码

=数据内容=

1. <xarray.Dataset>
   
2. Dimensions:    (longitude: 101, latitude: 61, level: 5, time: 144)
   
3. Coordinates:
   
4.   * longitude  (longitude) float32 100.0 100.2 100.5 100.8 ... 124.5 124.8 125.0
   
5.   * latitude   (latitude) float32 30.0 29.75 29.5 29.25 ... 15.5 15.25 15.0
   
6.   * level      (level) int32 500 700 850 925 1000
   
7.   * time       (time) datetime64[ns] 2023-09-07T08:00:00 ... 2023-09-13T07:00:00
   
8. Data variables:
   
9.     q          (time, level, latitude, longitude) float32 ...
   
10.     u          (time, level, latitude, longitude) float32 ...
    
11.     v          (time, level, latitude, longitude) float32 ...
    
12. Attributes:
    
13.     Conventions:  CF-1.6
    
14.     history:      2023-10-07 01:35:05 GMT by grib_to_netcdf-2.25.1: /opt/ecmw...

复制代码

1.2 数据处理与筛选*世界时UTC转北京时BJT

北京时比世界时快8个小时，选取数据内容的时间变量time。

1. def time_adjust(nc):

复制代码

=要素坐标=

1. data     = time_adjust(data)
   
2. time     = data.time
   
3. level    = data.level
   
4. lat      = data.latitude
   
5. lon      = data.longitude
   
6. lons,lats=np.meshgrid(lon,lat)
   
7. # time 选择第一个时间节点
   
8. time0   = time[0].values
   
9. # level 选择第四个高度节点
   
10. level0  = level[3].values

复制代码

=水汽及风速要素筛选=

1. u       = data.u.loc[time0,level0,:,:]  # 单位m/s
   
2. v       = data.v.loc[time0,level0,:,:]  # 单位m/s
   
3. q       = data.q.loc[time0,level0,:,:]  # 单位kg**-1

复制代码

二、绘图设置

1. import pandas as pd
   
2. import matplotlib.gridspec as gridspec
   
3. import metpy.constants as constants  # 里面是常数
   
4. &#8203;
   
5. import cartopy.crs as ccrs
   
6. import matplotlib.pyplot as plt
   
7. import cartopy.io.shapereader as shpreader
   
8. plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文
   
9. plt.rcParams['axes.unicode_minus']=False #用来正常显示负号

复制代码

=2.1 投影选择=

将地球三维球体投影到二维面上，减少失真。

主要方式有默认投影（PlateCarree）、兰勃脱投影（Lambert）、墨卡托投影（Mercator）、极投影。

参考 https://blog.csdn.net/qq_40497403/article/details/119248198

1. proj = ccrs.PlateCarree()
   
2. fig = plt.figure(figsize=(16, 12), dpi=100)
   
3. plotcrs = ccrs.PlateCarree()
   
4. ax = plt.subplot(111, projection=plotcrs)

复制代码

=2.2 设置范围=&#8203;

1. left   = 105
   
2. right  = 120
   
3. bottom = 16
   
4. upper  = 27
   
5. region = [left, right, bottom, upper]  # 要绘制的范围lon1,lon2,lat1,lat2
   
6. ax.set_extent(region, crs=proj)  # 设置绘图区范围！地球经纬度坐标

复制代码

*经纬度坐标设置

1. <blockquote>ax.set_xticks(np.arange(103, 120, 5), crs=ccrs.PlateCarree())  # 经度坐标轴

复制代码

=2.3 计算单层水汽通量和水汽通量散度=

g 的单位为m/s2，换算为N/kg，再换算为10-2hPa·m2/kg.

计算q*v/g,最终单层水汽通量的单位是10&#178;kg/m&#8226;hPa&#8226;s

10&#178;kg/m&#8226;hPa&#8226;s==10&#178;*1000g/100cm&#8226;hPa&#8226;s 转换成业务网的量纲单位 k g/cm&#8226;hPa&#8226;s

1. <blockquote># 计算方法

复制代码

三、要素分解
3.1 要素 ① 设置

采用的是LaTex数学公式，matplotlib是支持LaTex格式输入的。

参考：https://blog.csdn.net/fuchenqianxiao/article/details/124088512

1. # 右侧时间标签
   
2. time = pd.to_datetime(q.time.values)
   
3. hour = rename(time.hour)
   
4. day  = rename(time.day)
   
5. name = day + '日' + hour + '时' + rename(q.level.values) + 'hPa水汽通量图'
   
6. t3 = ax.set_title(r"{}/{}/{}/{}".format(time.year,time.month,time.day,rename(time.hour)),
   
7.                   loc='right', fontsize=15)
   
8. &#8203;
   
9. # 左侧备注标签
   
10. ## fontsize=20，总共20个字符，长度400
    
11. t1 = ax.set_title("Vapor Flux at {}hPa".format(level0),loc='left', fontsize=20)
    
12. t2 = ax.set_title(r'($\mathit{F}_{\mathit{H}}\mathit{:} \mathit{g/(cm·hPa·s)}$)', x=0.35,y=0.999, fontsize=12)3.2 要素 ② 设置
    
13. 色标上要素名称设置。
    
14. cb.ax.set_yticklabels([''] + [str(tick) for tick in cb.get_ticks()[1:-1]] + [''])
    
15. cb.ax.set_title(r'$\mathit{F}_{\mathit{H}}
    
16. [/align][align=left][b][b][size=5]3.3 要素 ③ 设置 [/size][/b][/b][align=left]见2.2经纬度坐标设置。[/align]
    
17. [b][size=4]*底图设置[/size][/b][align=left]底图的设置一定要在绘制填色图之后，否则会被填色图覆盖。[/align][code]china = shpreader.Reader(r'./basemap/bou_4l/bou2_4l.shp').geometries()
    
18. ax.add_geometries(china, plotcrs, facecolor='none',edgecolor='black', zorder=1)

复制代码

3.2 要素 ③ 设置

cb.ax.set_yticklabels([''] + [str(tick) for tick in cb.get_ticks()[1:-1]] + [''])

cb.ax.set_title(r'$\mathit{F}_{\mathit{H}}$', fontsize=15)

3.3 要素 ③ 设置

见2.2经纬度坐标设置。

*底图设置

底图的设置一定要在绘制填色图之后，否则会被填色图覆盖。

1. china = shpreader.Reader(r'./basemap/bou_4l/bou2_4l.shp').geometries()
   
2. ax.add_geometries(china, plotcrs, facecolor='none',edgecolor='black', zorder=1)

复制代码

3.4 要素 ④ 设置

绘制水汽通量风矢。

1. <blockquote># X,Y,U,V 确定位置和对应的风速

复制代码

四、完整代码

1. <blockquote>import xarray as xr

复制代码

沐雨淋觞

什么奇奇怪怪的排版工具，还找不到删帖在哪

沐雨淋觞

沐雨淋觞 发表于 2024-6-24 12:14
什么奇奇怪怪的排版工具，还找不到删帖在哪

算了，吞掉的代码不改了。https://blog.csdn.net/weixin_484 ... 1001.2014.3001.5501