WRF tools

Let’s open a WRF model output file with xarray first:

In [1]: import xarray as xr

In [2]: from salem.utils import get_demo_file

In [3]: ds = xr.open_dataset(get_demo_file('wrfout_d01.nc'))

WRF files are a bit messy:

In [4]: ds
Out[4]: 
<xarray.Dataset>
Dimensions:           (Time: 3, bottom_top: 27, bottom_top_stag: 28, south_north: 150, south_north_stag: 151, west_east: 150, west_east_stag: 151)
Coordinates:
    XLONG             (Time, south_north, west_east) float32 70.7231 70.9749 ...
    XLONG_U           (Time, south_north, west_east_stag) float32 70.5974 ...
    XLAT_U            (Time, south_north, west_east_stag) float32 7.76865 ...
    XLAT_V            (Time, south_north_stag, west_east) float32 7.66444 ...
    XLONG_V           (Time, south_north_stag, west_east) float32 70.7436 ...
    XLAT              (Time, south_north, west_east) float32 7.78907 7.82944 ...
  * Time              (Time) int64 0 1 2
  * south_north       (south_north) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 ...
  * west_east         (west_east) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 ...
  * bottom_top        (bottom_top) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 ...
  * west_east_stag    (west_east_stag) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 ...
  * south_north_stag  (south_north_stag) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 ...
  * bottom_top_stag   (bottom_top_stag) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 ...
Data variables:
    Times             (Time) |S19 '2008-10-26_12:00:00' ...
    T2                (Time, south_north, west_east) float32 301.977 302.008 ...
    RAINC             (Time, south_north, west_east) float32 0.0 0.0 0.0 0.0 ...
    RAINNC            (Time, south_north, west_east) float32 0.0 0.0 0.0 0.0 ...
    U                 (Time, bottom_top, south_north, west_east_stag) float32 1.59375 ...
    V                 (Time, bottom_top, south_north_stag, west_east) float32 -2.07031 ...
    PH                (Time, bottom_top_stag, south_north, west_east) float32 0.0 ...
    PHB               (Time, bottom_top_stag, south_north, west_east) float32 0.0 ...
Attributes:
    note: Global attrs removed.

WRF files aren’t exactly CF compliant: you’ll need a special parser for the timestamp, the coordinate names are a bit exotic and do not correspond to the dimension names, they contain so-called staggered variables (and their correponding coordinates), etc.

Salem defines a special parser for these files:

In [5]: import salem

In [6]: ds = salem.open_wrf_dataset(get_demo_file('wrfout_d01.nc'))

This parser greatly simplifies the file structure:

In [7]: ds
Out[7]: 
<xarray.Dataset>
Dimensions:       (bottom_top: 27, south_north: 150, time: 3, west_east: 150)
Coordinates:
    lon           (south_north, west_east) float32 70.7231 70.9749 71.2267 ...
    lat           (south_north, west_east) float32 7.78907 7.82944 7.86932 ...
  * time          (time) datetime64[ns] 2008-10-26T12:00:00 ...
  * south_north   (south_north) float64 -2.235e+06 -2.205e+06 -2.175e+06 ...
  * west_east     (west_east) float64 -2.235e+06 -2.205e+06 -2.175e+06 ...
  * bottom_top    (bottom_top) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ...
Data variables:
    T2            (time, south_north, west_east) float32 301.977 302.008 ...
    RAINC         (time, south_north, west_east) float32 0.0 0.0 0.0 0.0 0.0 ...
    RAINNC        (time, south_north, west_east) float32 0.0 0.0 0.0 0.0 0.0 ...
    U             (time, bottom_top, south_north, west_east) float32 1.73438 ...
    V             (time, bottom_top, south_north, west_east) float32 -2.125 ...
    PH            (time, bottom_top, south_north, west_east) float32 18.7031 ...
    PHB           (time, bottom_top, south_north, west_east) float32 277.504 ...
    T2C           (time, south_north, west_east) float32 28.8266 28.8578 ...
    PRCP          (time, south_north, west_east) float32 nan nan nan nan nan ...
    WS            (time, bottom_top, south_north, west_east) float32 2.74293 ...
    GEOPOTENTIAL  (time, bottom_top, south_north, west_east) float32 296.207 ...
    Z             (time, bottom_top, south_north, west_east) float32 30.1944 ...
    PRCP_NC       (time, south_north, west_east) float32 nan nan nan nan nan ...
    PRCP_C        (time, south_north, west_east) float32 nan nan nan nan nan ...
Attributes:
    note: Global attrs removed.
    pyproj_srs: +units=m +proj=lcc +lat_1=29.0400009155 +lat_2=29.0400009155 +lat_0=29.0399971008 +lon_0=89.8000030518 +x_0=0 +y_0=0 +a=6370000 +b=6370000

Note that some dimensions / coordinates have been renamed, new variables have been defined, and the staggered dimensions have disappeared.

Diagnostic variables

Salem adds a layer between xarray and the underlying NetCDF file. This layer computes new variables “on the fly” or, in the case of staggered variables, “unstaggers” them:

In [8]: ds.U
Out[8]: 
<xarray.DataArray 'U' (time: 3, bottom_top: 27, south_north: 150, west_east: 150)>
[1822500 values with dtype=float32]
Coordinates:
    lon          (south_north, west_east) float32 70.7231 70.9749 71.2267 ...
    lat          (south_north, west_east) float32 7.78907 7.82944 7.86932 ...
  * time         (time) datetime64[ns] 2008-10-26T12:00:00 ...
  * south_north  (south_north) float64 -2.235e+06 -2.205e+06 -2.175e+06 ...
  * west_east    (west_east) float64 -2.235e+06 -2.205e+06 -2.175e+06 ...
  * bottom_top   (bottom_top) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ...
Attributes:
    FieldType: 104
    MemoryOrder: XYZ
    description: x-wind component
    units: m s-1
    stagger: X
    pyproj_srs: +units=m +proj=lcc +lat_1=29.0400009155 +lat_2=29.0400009155 +lat_0=29.0399971008 +lon_0=89.8000030518 +x_0=0 +y_0=0 +a=6370000 +b=6370000

This computation is done only on demand (just like a normal NetCDF variable), this new layer is therefore relatively cheap.

In addition to unstaggering, Salem adds a number of “diagnostic” variables to the dataset. Some are convenient (like T2C, temperature in Celsius instead of Kelvins), but others are more complex (e.g. SLP for sea-level pressure, or PRCP which computes step-wize total precipitation out of the accumulated fields). For a list of diagnostic variables (and TODOs!), refer to GH18.

In [9]: ds.PRCP.isel(time=-1).salem.quick_map(cmap='Purples', vmax=5)
Out[9]: <salem.graphics.Map at 0x7fd3c0ad6910>
_images/plot_wrf_diag.png

Vertical interpolation

The WRF vertical coordinates are eta-levels, which is not a very practical coordinate for analysis or plotting. With the functions wrf_zlevel() and wrf_plevel() it is possible to interpolate the 3d data at either altitude or pressure levels:

In [10]: ws_h = ds.isel(time=1).salem.wrf_zlevel('WS', levels=10000.)

In [11]: ws_h.salem.quick_map();
_images/plot_wrf_zinterp.png

Note that currently, the interpolation is quite slow, see GH25. It’s probably wise to do it on single time slices or aggregated data, rather than huge data cubes.

Geogrid simulator

Salem provides a small tool which comes handy when defining new WRF domains. It parses the WPS configuration file (namelist.wps) and generates the grids and maps corresponding to each domain.

geogrid_simulator() will search for the &geogrid section of the file and parse it:

In [12]: fpath = get_demo_file('namelist_mercator.wps')

In [13]: with open(fpath, 'r') as f:  # this is just to show the file
   ....:     print(f.read())
   ....: 
&geogrid
 parent_id         =   1, 1, 2, 3,
 parent_grid_ratio =   1, 3, 5, 5,
 i_parent_start    =   1, 28, 54, 40,
 j_parent_start    =   1, 24, 97, 50,
 e_we              =  211, 178, 101, 136,
 e_sn              =  131, 220, 141, 136,
 dx = 30000,
 dy = 30000,
 map_proj = 'mercator',
 ref_lat   =  -8.0,
 ref_lon   =  -62.0,
 truelat1  =  -8.0,
/



In [14]: from salem import geogrid_simulator

In [15]: g, maps = geogrid_simulator(fpath)

In [16]: maps[0].set_rgb(natural_earth='lr')  # add a background image

In [17]: maps[0].visualize(title='Domains 1 to 4')
_images/plot_geo_simu.png