Reading WRF model data

Basemap is specially good at drawing numerical weather prediction models outputs, such as WRF. WRF is a widely used model, but most of the example can be run with other model outputs, just adapting the variable names.

At the UCAR website is possible to download a WRF sample output file

The output file descriptor (cdl) contains all the information about the model size, fields, projection, etc. For instance, we will need the projection information to project the output properly:

:CEN_LAT = 34.83158f ;

:CEN_LON = -81.02756f ;

:TRUELAT1 = 30.f ;

:TRUELAT2 = 60.f ;

Note

Make sure that the gdal library installed can read NetCDF files. Check it writing gdalinfo –formats

Plotting a field as a contour

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()
../_images/read_wrf.png

  • Note how does GDAL open the files. When dealing with NetCDF files, it uses what it calls subdatasets, so every variable acts as an independent file

  • XLONG and XLAT contain the information about the longitude and latitude of every point in the matrix. This is very useful using Basemap, since the fields needed for the methods to indicate those positions are already calculated.

  • When drawing the contour, the first band from each variable is taken, since in this case, the model has several bands for longitudes, latitudes and temperature

  • The strange shape of the data is because the output map is in the default Basemap projection (Equirectangular), but the model is not, so it has to be re-projected.

Projecting the map

As its name suggests, the model uses a Lambert conformal projection. The metadata file tells us how is this Lambert Conformal projection defined, so plotting the map is easy:

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()
../_images/read_wrf_projected.png

  • The limits in the map don’t match those of the data to show that now, the shape of the data is rectangular, so the projection is properly configured

Wind barbs

Basemap has a function that makes drawing wind barbs very simple (unlike using other GIS libraries, btw):

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()
../_images/read_wrf_barbs.png

  • Now, the limits of the map match those of the data, so no white space is shown. The resolution parameter is also changed, to get a prettier coastline.

  • Not all the possible barbs are drawn, since the map would be difficult to understand.

    • To eliminate some of the points, the numpy.arange function is used, to select the pixels to be used

    • After having to arrays with the positions, a 2d matrix is created with this values, using numpy.meshgrid. Now is possible to select only these points from any array, as shown when using the barbs method

  • A contour with the wind speed is plotted under the barbs, to make the map more interesting. Note that since the fields are numpy arrays, the module is easily calculated