Basemap系列教程：读取WRF模式数据

bugsuse

发布于 2020-04-20 13:47:33

1.8K1

发布于 2020-04-20 13:47:33

文章被收录于专栏：气象杂货铺气象杂货铺

Basemap 特别擅长绘制数值天气模式输出数据，比如 WRF。WRF [注1] 模式是广泛使用的数值预报模式，只要变量名合适，大部分情况下都可以使用其它模式的输出来运行。

在 UCAR 网站可以下载 WRF 输出数据 [注2]。

输出文件描述 [注3] 包含了模式规模，域，投影等信息。例如，我们需要投影相关的信息来正确投影输出。

:CEN_LAT = 34.83158f ;

:CEN_LON = -81.02756f ;

:TRUELAT1 = 30.f ;

:TRUELAT2 = 60.f ;

注意：

确保 gdal 库安装了，从而便于读取 NetCDF 文件。

译注：读取 NetCDF文件也可以使用 netcdf4 库

绘制域

from osgeo import gdal
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

wrf_out_file = "wrfout_v2_Lambert.nc"

ds_lon = gdal.Open('NETCDF:"'+wrf_out_file+'":XLONG')
ds_lat = gdal.Open('NETCDF:"'+wrf_out_file+'":XLAT')
ds_t2 = gdal.Open('NETCDF:"'+wrf_out_file+'":T2')

map = Basemap(llcrnrlon=-95.,llcrnrlat=24.,urcrnrlon=-66.,urcrnrlat=45.)

map.contourf(ds_lon.ReadAsArray()[1], ds_lat.ReadAsArray()[1], ds_t2.ReadAsArray()[1]) 

map.drawcoastlines()
 
plt.show()

注意 GDAL 是如何读取 NetCDF文件的。读取文件时，它使用的是调用 subdatasets 即，子数据集的方式，因此每个变量都像是单独的文件
XLONG 和 XLAT 包含的是矩阵中每个点的经纬度信息。使用 Basemap 时这些信息是非常重要的
绘图时，仅使用每个变量的第一条带，这是因为longitudes， latitudes 和temperature 都有多条带（译注：因为每个变量都有多个时刻的输出）
由于 Basemap 默认的投影和模式输出采用的投影方式不同，导致结果看起来很奇怪。因此需要重新投影

投影地图

模式使用的是 Lambert conformal 投影，因此 Basemap 应该也设置为 Lambert conformal 。

from osgeo import gdal
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

wrf_out_file = "wrfout_v2_Lambert.nc"

ds_lon = gdal.Open('NETCDF:"'+wrf_out_file+'":XLONG')
ds_lat = gdal.Open('NETCDF:"'+wrf_out_file+'":XLAT')

ds_t2 = gdal.Open('NETCDF:"'+wrf_out_file+'":T2')

map = Basemap(llcrnrlon=-95.,llcrnrlat=27.,urcrnrlon=-65.,urcrnrlat=40.,
              projection='lcc', lat_1=30.,lat_2=60.,lat_0=34.83158,lon_0=-98.)

x, y = map(ds_lon.ReadAsArray()[1], ds_lat.ReadAsArray()[1])

map.contourf(x, y, ds_t2.ReadAsArray()[1]) 

map.drawcoastlines()
plt.show()

地图的边界和数据的边界并没有很好的匹配

Wind barbs

Basemap 有函数可以很容易绘制 wind barbs（不像其它 gis 库）

from osgeo import gdal
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np

wrf_out_file = "wrfout_v2_Lambert.nc"

ds_lon = gdal.Open('NETCDF:"'+wrf_out_file+'":XLONG')
ds_lat = gdal.Open('NETCDF:"'+wrf_out_file+'":XLAT')

ds_u = gdal.Open('NETCDF:"'+wrf_out_file+'":U10')
ds_v = gdal.Open('NETCDF:"'+wrf_out_file+'":V10')

map = Basemap(llcrnrlon=-93.7, llcrnrlat=28., urcrnrlon=-66.1, urcrnrlat=39.5,
              resolution = 'l',
              projection='lcc', lat_1=30., lat_2=60., lat_0=34.83158, lon_0=-98.)

x, y = map(ds_lon.ReadAsArray()[1], ds_lat.ReadAsArray()[1])


u = ds_u.ReadAsArray()[1]
v = ds_v.ReadAsArray()[1]

yy = np.arange(0, y.shape[0], 4)
xx = np.arange(0, x.shape[1], 4)

points = np.meshgrid(yy, xx) 

map.contourf(x, y, np.sqrt(u*u + v*v), alpha = 0.4) 
map.barbs(x[points], y[points], u[points], v[points]) 

map.drawcoastlines(color = '0.15')
plt.show()

现在，地图的边界和数据边界一致了，也没有空白区域了。resolution 参数也改变了，使海岸线更精细了
没有绘制所有的 barbs ，否则地图看起来会很难理解 1）使用 numpy.arange 函数过滤一些点 2）使用 meshgrid 函数创建经纬度坐标，后面进行绘图
barbs 图上叠加了风速等值线图，使图看起来更有意思。注意：由于这些变量都是 numpy 数组，因此计算起来很容易。

注1：

https://en.wikipedia.org/wiki/Weather_Research_and_Forecasting_Model

注2：

http://www.unidata.ucar.edu/software/netcdf/examples/wrfout_v2_Lambert.nc

注3：

https://www.unidata.ucar.edu/software/netcdf/examples/wrfout_v2_Lambert.cdl

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2017-03-17，如有侵权请联系 cloudcommunity@tencent.com 删除

编程算法