Python气象数据处理与绘图：基于数值模式资料和站点观测数据绘制相对湿度分布图

气海无涯 2021-07-12

5172

目的和意义：如何引入网格化数据之外的其他信息，将站点数据和网格化数据结合在一起，获取更高精度、更高分辨率和更高准确度的气象要素分布图。

实用性：格点化的天气预报数据往往数个小时前就已经获得，而再分析资料现在暂未获得。我们有的只有十几分钟前获得的站点观测数据，而我们需要的是一个较为准确的区域的当前网格化观测数据。

主要步骤：

下载GFS数值模式的预报数据和地面自动气象站观测数据
读取GFS预报数据的地表温度和露点温度变量
基于地表温度和露点温度，计算网格化地表相对湿度
使用cressman插值，以站点经纬度为第一猜测场，引入站点相对湿度，获得一个新的相对湿度
将GFS预报的相对湿度由经纬度网格处理为WRF的兰伯特投影

需要提前安装模块：siphon、metpy、cmaps、cartopy

一、下载GFS预报数据

1from siphon.catalog import TDSCatalog
2best_gfs = TDSCatalog('http://thredds.ucar.edu/thredds/catalog/grib/NCEP/GFS/'
3                      'Global_0p25deg/catalog.xml?dataset=grib/NCEP/GFS/Global_0p25deg/Best')
4best_ds = list(best_gfs.datasets.values())[0]
5ncss = best_ds.subset()
6query = ncss.query()

1from datetime import datetime, timedelta
2mytime = datetime.strptime('2020-08-13T00:00:00Z','%Y-%m-%dT%H:%M:%Sz')
3query.lonlat_box(north=55, south=0, east=140, west=70).time(mytime)
4query.variables('Temperature_surface','Dewpoint_temperature_height_above_ground',
5                'Relative_humidity_height_above_ground')
6query.accept('netcdf4')

输出：

1>>>var=Relative_humidity_height_above_ground&var=Temperature_surface&var=Dewpoint_temperature_height_above_ground&time=2020-08-13T00%3A00%3A00&west=70&east=140&south=0&north=55&accept=netcdf4

1from xarray.backends import NetCDF4DataStore
2import xarray as xr
3nc = ncss.get_data(query)
4data = xr.open_dataset(NetCDF4DataStore(nc))

1data

输出：

1data.to_netcdf('GFS_2020081300.nc','w')

1import shutil
2shutil.move('GFS_2020081300.nc', './work')

输出：

1'./work/GFS_2020081300.nc'

注释：以上仅介绍数据下载的方法，后续代码所用数据为2020年8月10日00时(UTC)的GFS数据。

二、读取GFS预报数据

注释：这里要注意维度问题，因为需要的是某一时刻的数据，所以后面把时间维度删去了。

1import metpy
2import xarray as xr
3data = xr.open_dataset('./work/GFS_2020081000.nc')
4data

1tem_2m = data['Temperature_surface'].metpy.unit_array.squeeze()
2dewpoint_2m = data['Dewpoint_temperature_height_above_ground'].metpy.unit_array.squeeze()
3rh_2m = data['Relative_humidity_height_above_ground'].metpy.unit_array.squeeze()
4lon_1d = data['lon']
5lat_1d = data['lat']
6data.close()

三、计算相对湿度(基于地表温度和露点温度)

注释：GFS预报场其实是包含相对湿度变量的，根据作业要求，这里基于地表温度和露点温度，计算网格化地表相对湿度。在后面也会比较计算结果和预报结果的差异。计算参考：相对湿度计算说明。发现算出来的相对湿度有大于1的情况，这是由于水汽过饱和导致的，后面将大于1的全部当做1处理。

1import metpy.calc
2rh_2m_cal = metpy.calc.relative_humidity_from_dewpoint(tem_2m, dewpoint_2m)

1print(np.nanmax(rh_2m_cal))
2>>>1.70308518409729 dimensionless

1print(np.nanmin(rh_2m_cal))
2>>>0.09976071119308472 dimensionless

1rh_2m_cal[rh_2m_cal>1]=1

1print(np.nanmax(rh_2m_cal))
2>>>1.0 dimensionless

1rh_2m_cal = rh_2m_cal*100

四、绘图展示数据

注释：由于数据原本是等经纬的，所以目前先不进行投影，仍保留原始状态

1lon_2d, lat_2d = np.meshgrid(lon_1d, lat_1d)

1%matplotlib inline
2import matplotlib.pyplot as plt
3import cartopy.crs as ccrs
4import cmaps
5import numpy as np
6import cartopy.io.shapereader as shpreader
7from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

 1#第一张图：GFS的2m气温预报场
 2fig = plt.figure(figsize=(15, 12))
 3ax = plt.axes(projection=ccrs.PlateCarree())
 4ax.set_extent([70, 140, 0, 50])
 5contours = ax.contourf(lon_2d, lat_2d, tem_2m.to('degC'), levels=range(-5,46), 
 6                       transform=ccrs.PlateCarree(),cmap=cmaps.NCV_bright)
 7province = shpreader.Reader('./work/province.shp')
 8nineline = shpreader.Reader('./work/nine_line.shp')
 9ax.add_geometries(province.geometries(), crs=ccrs.PlateCarree(), linewidths=0.5,edgecolor='k',facecolor='none')
10ax.add_geometries(nineline.geometries(), crs=ccrs.PlateCarree(),linewidths=0.5, color='k')
11gl = ax.gridlines(draw_labels=True, linewidth=1, color='k', alpha=0.5, linestyle='--')
12gl.xlabels_top = False  #关闭顶端标签
13gl.ylabels_right = False  #关闭右侧标签
14gl.xformatter = LONGITUDE_FORMATTER  #x轴设为经度格式
15gl.yformatter = LATITUDE_FORMATTER  #y轴设为纬度格式
16cbar = fig.colorbar(contours,shrink=0.8)
17cbar.set_label('Temperature(℃)',fontsize=15)  
18cbar.set_ticks(np.arange(-5,46,5))
19ax.set_title('GFS_Temperature_2m  valid at 2020-08-10T00:00:00(UTC)',fontsize=15)
20plt.show()

 1#第二张图：GFS的2m露点预报场
 2fig = plt.figure(figsize=(15, 12))
 3ax = plt.axes(projection=ccrs.PlateCarree())
 4ax.set_extent([70, 140, 0, 50])
 5contours = ax.contourf(lon_2d, lat_2d, dewpoint_2m.to('degC'), levels=range(-15,36), 
 6                       transform=ccrs.PlateCarree(),cmap=cmaps.NCV_bright)
 7province = shpreader.Reader('./work/province.shp')
 8nineline = shpreader.Reader('./work/nine_line.shp')
 9ax.add_geometries(province.geometries(), crs=ccrs.PlateCarree(), linewidths=0.5,edgecolor='k',facecolor='none')
10ax.add_geometries(nineline.geometries(), crs=ccrs.PlateCarree(),linewidths=0.5, color='k')
11gl = ax.gridlines(draw_labels=True, linewidth=1, color='k', alpha=0.5, linestyle='--')
12gl.xlabels_top = False  #关闭顶端标签
13gl.ylabels_right = False  #关闭右侧标签
14gl.xformatter = LONGITUDE_FORMATTER  #x轴设为经度格式
15gl.yformatter = LATITUDE_FORMATTER  #y轴设为纬度格式
16cbar = fig.colorbar(contours,shrink=0.8)
17cbar.set_label('Dewpoint Temperature(℃)',fontsize=15)  
18cbar.set_ticks(np.arange(-15,36,5))
19ax.set_title('GFS_Dewpoint_Temperature_2m  valid at 2020-08-10T00:00:00(UTC)',fontsize=15)
20plt.show()

 1#第三张图：GFS的2m相对湿度预报场
 2fig = plt.figure(figsize=(15, 12))
 3ax = plt.axes(projection=ccrs.PlateCarree())
 4ax.set_extent([70, 140, 0, 50])
 5contours = ax.contourf(lon_2d, lat_2d, rh_2m, levels=range(0,101), 
 6                       transform=ccrs.PlateCarree(),cmap=cmaps.MPL_terrain_r)
 7province = shpreader.Reader('./work/province.shp')
 8nineline = shpreader.Reader('./work/nine_line.shp')
 9ax.add_geometries(province.geometries(), crs=ccrs.PlateCarree(), linewidths=0.5,edgecolor='k',facecolor='none')
10ax.add_geometries(nineline.geometries(), crs=ccrs.PlateCarree(),linewidths=0.5, color='k')
11gl = ax.gridlines(draw_labels=True, linewidth=1, color='k', alpha=0.5, linestyle='--')
12gl.xlabels_top = False  #关闭顶端标签
13gl.ylabels_right = False  #关闭右侧标签
14gl.xformatter = LONGITUDE_FORMATTER  #x轴设为经度格式
15gl.yformatter = LATITUDE_FORMATTER  #y轴设为纬度格式
16cbar = fig.colorbar(contours,shrink=0.8)
17cbar.set_label('Relative Humidity(%)',fontsize=15)  
18cbar.set_ticks(np.arange(0,101,5))
19ax.set_title('GFS_Relative Humidity_2m  valid at 2020-08-10T00:00:00(UTC)',fontsize=15)
20plt.show()

 1#第四张图：GFS的2m相对湿度预报场（计算所得）
 2fig = plt.figure(figsize=(15, 12))
 3ax = plt.axes(projection=ccrs.PlateCarree())
 4ax.set_extent([70, 140, 0, 50])
 5contours = ax.contourf(lon_2d, lat_2d, rh_2m_cal, levels=range(0,101), 
 6                       transform=ccrs.PlateCarree(),cmap=cmaps.MPL_terrain_r)
 7province = shpreader.Reader('./work/province.shp')
 8nineline = shpreader.Reader('./work/nine_line.shp')
 9ax.add_geometries(province.geometries(), crs=ccrs.PlateCarree(), linewidths=0.5,edgecolor='k',facecolor='none')
10ax.add_geometries(nineline.geometries(), crs=ccrs.PlateCarree(),linewidths=0.5, color='k')
11gl = ax.gridlines(draw_labels=True, linewidth=1, color='k', alpha=0.5, linestyle='--')
12gl.xlabels_top = False  #关闭顶端标签
13gl.ylabels_right = False  #关闭右侧标签
14gl.xformatter = LONGITUDE_FORMATTER  #x轴设为经度格式
15gl.yformatter = LATITUDE_FORMATTER  #y轴设为纬度格式
16cbar = fig.colorbar(contours,shrink=0.8)
17cbar.set_label('Relative Humidity(%)',fontsize=15)  
18cbar.set_ticks(np.arange(0,101,5))
19ax.set_title('GFS_Relative Humidity_2m[calculated]  valid at 2020-08-10T00:00:00(UTC)',fontsize=15)
20plt.show()

五、自动气象站观测数据(站点资料)

注释：此处数据来自中国气象数据网，这里展示的数据并没有包含全国所有自动站，华南和东南地区没有下载，其实可以看看站点空间分布不均匀带来的影响。

1import pandas as pd
2
3df = pd.read_excel('./work/20200810_RH.xls')
4df

1station_lon = df['经度']
2station_lat = df['纬度']
3station_RH = df['相对湿度']

 1#第五张图：自动站的相对湿度观测资料(散点图)
 2fig = plt.figure(figsize=(15, 12))
 3ax = plt.axes(projection=ccrs.PlateCarree())
 4ax.set_extent([70, 140, 0, 50])
 5sc = ax.scatter(station_lon, station_lat, c=station_RH, s=5, marker='o', vmin=0, vmax=100, 
 6                cmap=cmaps.GMT_drywet,transform=ccrs.PlateCarree())
 7province = shpreader.Reader('./work/province.shp')
 8nineline = shpreader.Reader('./work/nine_line.shp')
 9ax.add_geometries(province.geometries(), crs=ccrs.PlateCarree(), linewidths=0.5,edgecolor='k',facecolor='none')
10ax.add_geometries(nineline.geometries(), crs=ccrs.PlateCarree(),linewidths=0.5, color='k')
11gl = ax.gridlines(draw_labels=True, linewidth=1, color='k', alpha=0.5, linestyle='--')
12gl.xlabels_top = False  #关闭顶端标签
13gl.ylabels_right = False  #关闭右侧标签
14gl.xformatter = LONGITUDE_FORMATTER  #x轴设为经度格式
15gl.yformatter = LATITUDE_FORMATTER  #y轴设为纬度格式
16ax.stock_img()
17cbar = fig.colorbar(sc,shrink=0.8)
18cbar.set_label('Relative Humidity(%)',fontsize=15)  
19cbar.set_ticks(np.arange(0,101,5))
20ax.set_title('OBS_Relative Humidity  valid at 2020-08-10T00:00:00(UTC)',fontsize=15)
21plt.show()

六、定义cressman插值函数并进行插值

注释：这部分是最核心的部分，以下代码来自冬青老师的课件，太有用了。因为计算量较大，使用了numba加速，不然效率很低。但后面找到了一个更好的办法。

 1from math import acos, asin, cos, pi, radians, sin, sqrt
 2from numba import *
 3import multiprocessing
 4import warnings
 5warnings.filterwarnings('ignore')
 6#定义cressman插值函数
 7def cressman_interpolation(station_lat,station_lon,station_data,lat_line,lon_line, search_radius,first_guess=None):
 8    scale = 1.  # 校正系数
 9    grid_i ,grib_j = lon_grid.shape
10    if first_guess is None:
11        first_guess = np.zeros([grid_i ,grib_j], np.float32) 
12    result = np.zeros([grid_i ,grib_j], np.float32)    #最终结果值
13    station_number = len(station_lat)  #站点个数 
14    # 插值，计算网格在每个观测站点上的值
15    station_interp_data = np.zeros((station_number))
16    first_guess_xr = xr.DataArray(first_guess, coords=[lat_line, lon_line], dims=[ 'lat','lon'])
17    for istation in range(station_number):
18        station_interp_data[istation] = first_guess_xr.interp(lat=station_lat[istation],lon= station_lon[istation]).values  #站点的第一猜测值
19        # print(station_interp_data[istation])
20    print("station_interp_data generated")
21    ################################
22    lat_number,lon_number = lon_grid.shape
23    result = cressman_gradually_correction(lat_number,lon_number,station_number,station_lat,station_lon,lat_grid,lon_grid,search_radius,station_data,station_interp_data,result)
24    return result
25@jit
26def cressman_gradually_correction(lat_number,lon_number,station_number,station_lat,station_lon,lat_grid,lon_grid,search_radius,station_data,station_interp_data,result):
27    for ilat in range(lat_number):
28        for ilon in range(lon_number):
29            sum_w2deltaa = 0.0
30            sum_w = 0.0
31            for istation in range(station_number):
32                #计算每个台站和格点的距离
33                # 先计算大致距离（曼哈顿距离），减少计算量
34                manhattan_distance = abs(station_lat[istation]-lat_grid[ilat,ilon])+abs(station_lon[istation]-lon_grid[ilat,ilon])
35                if manhattan_distance*111 > search_radius*2/1.3:  #
36                    w = 0  #权重系数
37                else:
38                    distance = geodistance(station_lat[istation],station_lon[istation],lat_grid[ilat,ilon],lon_grid[ilat,ilon])
39                    if distance >= search_radius:
40                        w = 0  #权重系数
41                    else:
42                        w = (search_radius ** 2 - distance ** 2) / (search_radius ** 2 + distance ** 2)
43                        sum_w += w 
44                        sum_w2deltaa += w**2*(station_data[istation]-station_interp_data[istation])
45
46            if sum_w > 0:
47                # print("lat,lon=",ilat,ilon)
48                # print("sum_w=",sum_w)
49                deltaa = sum_w2deltaa / sum_w    
50                # print("deltaa=",deltaa)
51                result[ilat,ilon] +=  deltaa
52    return result
53@jit
54def geodistance(lat1,lng1,lat2,lng2):  #计算地球上两点距离
55    lng1, lat1, lng2, lat2 = map(radians, [float(lng1), float(lat1), float(lng2), float(lat2)]) # 经纬度转换成弧度
56    dlon=lng2-lng1
57    dlat=lat2-lat1
58    a=sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
59    distance=2*asin(sqrt(a))*6371 # 地球平均半径，6371km
60    return distance

1#cressman插值
2from metpy.interpolate import interpolate_to_grid, remove_nan_observations, remove_repeat_coordinates
3station_lon,station_lat,station_RH = remove_nan_observations(station_lon,station_lat,station_RH)
4station_lon,station_lat,station_RH = remove_repeat_coordinates(station_lon,station_lat,station_RH)
5lon_grid, lat_grid, first_guess = interpolate_to_grid(station_lon,station_lat,station_RH, interp_type='cressman', minimum_neighbors=1,
6                                  hres=0.25, search_radius=0.5)

 1#第六张图：自动站的相对湿度空间插值
 2fig = plt.figure(figsize=(15, 12))
 3ax = plt.axes(projection=ccrs.PlateCarree())
 4ax.set_extent([70, 140, 0, 50])
 5contours = ax.contourf(lon_grid, lat_grid, first_guess, levels=range(0,101), 
 6                       transform=ccrs.PlateCarree(),cmap=cmaps.MPL_terrain_r)
 7province = shpreader.Reader('./work/province.shp')
 8nineline = shpreader.Reader('./work/nine_line.shp')
 9ax.add_geometries(province.geometries(), crs=ccrs.PlateCarree(), linewidths=0.5,edgecolor='k',facecolor='none')
10ax.add_geometries(nineline.geometries(), crs=ccrs.PlateCarree(),linewidths=0.5, color='k')
11gl = ax.gridlines(draw_labels=True, linewidth=1, color='k', alpha=0.5, linestyle='--')
12gl.xlabels_top = False  #关闭顶端标签
13gl.ylabels_right = False  #关闭右侧标签
14gl.xformatter = LONGITUDE_FORMATTER  #x轴设为经度格式
15gl.yformatter = LATITUDE_FORMATTER  #y轴设为纬度格式
16cbar = fig.colorbar(contours,shrink=0.8)
17cbar.set_label('Relative Humidity(%)',fontsize=15)  
18cbar.set_ticks(np.arange(0,101,5))
19ax.set_title('Guess_Relative Humidity_2m  valid at 2020-08-10T00:00:00(UTC)',fontsize=15)
20plt.show()

注释：分辨率和搜索半径很重要，这个结果比较可靠，因为站点所能代表的范围有限，当然还应该考虑地形、土地利用等信息才能得到更为准确的空间分布。

七、简单的数据同化(数据融合)

目前已经得到两个相对湿度的格点数据了，分别是GFS预报场和地面站点的观测场，问题是预报场可靠性存疑，但空间上连续，而观测场可靠性较高，但空间上有缺失。这里需要将两者结合起来，得到既可靠又空间连续的相对湿度分布。空间分辨率：（1）GFS预报场：0.25×0.25；（2）地面观测场：0.25×0.25。如果遇到非一致的空间分辨率还得进行插值，使两个空间属性一致。但是，除了空间分辨率问题，还有格点的一一对应关系，所以这里还是要插值qaq。

1first_guess.shape
2
3>>>(126, 237)
4
5rh_2m_cal.shape
6
7>>>(221, 281)

1lon_2d, lat_2d = np.meshgrid(lon_1d, lat_1d)
2lon_2d = xr.DataArray(lon_2d, dims=['x','y'])
3lat_2d = xr.DataArray(lat_2d, dims=['x','y'])

1lon_grid = xr.DataArray(lon_grid, dims=['x','y'])
2lat_grid = xr.DataArray(lat_grid, dims=['x','y'])

1first_guess_xr= xr.DataArray(first_guess, coords=[lat_grid[:,0], lon_grid[0]], dims=["lat", "lon"])

1first_guess_xr

1first_guess_new = first_guess_xr.interp(lon=lon_2d, lat=lat_2d)
2
3first_guess_new

 1#第七张图：重新插值的first_guess
 2fig = plt.figure(figsize=(15, 12))
 3ax = plt.axes(projection=ccrs.PlateCarree())
 4ax.set_extent([70, 140, 0, 50])
 5contours = ax.contourf(lon_2d, lat_2d, first_guess_new, levels=range(0,101), 
 6                       transform=ccrs.PlateCarree(),cmap=cmaps.MPL_terrain_r)
 7province = shpreader.Reader('./work/province.shp')
 8nineline = shpreader.Reader('./work/nine_line.shp')
 9ax.add_geometries(province.geometries(), crs=ccrs.PlateCarree(), linewidths=0.5,edgecolor='k',facecolor='none')
10ax.add_geometries(nineline.geometries(), crs=ccrs.PlateCarree(),linewidths=0.5, color='k')
11gl = ax.gridlines(draw_labels=True, linewidth=1, color='k', alpha=0.5, linestyle='--')
12gl.xlabels_top = False  #关闭顶端标签
13gl.ylabels_right = False  #关闭右侧标签
14gl.xformatter = LONGITUDE_FORMATTER  #x轴设为经度格式
15gl.yformatter = LATITUDE_FORMATTER  #y轴设为纬度格式
16cbar = fig.colorbar(contours,shrink=0.8)
17cbar.set_label('Relative Humidity(%)',fontsize=15)  
18cbar.set_ticks(np.arange(0,101,5))
19ax.set_title('first_guess_new  valid at 2020-08-10T00:00:00(UTC)',fontsize=15)
20plt.show()

1first_guess_new = np.array(first_guess_new)
2rh_2m_cal = np.array(rh_2m_cal)
3rh_2m_cal[~np.isnan(first_guess_new)] = first_guess_new[~np.isnan(first_guess_new)]

八、变为WRF的兰伯特投影(重投影)

 1#第八张图：GFS的2m相对湿度预报场（计算所得）
 2proj = ccrs.LambertConformal(central_longitude=105,central_latitude=90,standard_parallels=(25, 47))
 3fig = plt.figure(figsize=(15, 12))
 4ax = fig.subplots(1, 1, subplot_kw={'projection': proj}) 
 5ax.set_extent([70, 140, 0, 50])
 6contours = ax.contourf(lon_2d, lat_2d, rh_2m_cal, levels=range(0,101), 
 7                       transform=ccrs.PlateCarree(),cmap=cmaps.MPL_terrain_r)
 8province = shpreader.Reader('./work/province.shp')
 9nineline = shpreader.Reader('./work/nine_line.shp')
10ax.add_geometries(province.geometries(), crs=ccrs.PlateCarree(), linewidths=0.5,edgecolor='k',facecolor='none')
11ax.add_geometries(nineline.geometries(), crs=ccrs.PlateCarree(),linewidths=0.5, color='k')
12cbar = fig.colorbar(contours,shrink=0.8)
13cbar.set_label('Relative Humidity(%)',fontsize=15)  
14cbar.set_ticks(np.arange(0,101,5))
15ax.set_title('GFS_Relative Humidity_2m[assimilation]  valid at 2020-08-10T00:00:00(UTC)',fontsize=15)
16plt.show()

注释：这次发现经纬度标不上了，因为cartopy版本是0.17.0，大家升级到0.18.0就可以对兰伯特标注经纬度了。

有问题可以到QQ群里进行讨论，我们在那边等大家。

QQ群号：854684131

python

文章转载自气海无涯，如果涉嫌侵权，请发送邮件至：contact@modb.pro进行举报，并提供相关证据，一经查实，墨天轮将立刻删除相关内容。