wrf-python/test/ipynb/WRF_Workshop_Demo.ipynb

Demo for wrf-python¶

1.0 Introduction¶

1) Created to provide Python with the same functionality found in WRF-NCL

* WRF variable computation and extraction
* 2D and 3D interpolation routines
* WRF-ARW only
* Metadata

2) Key differences with WRF-NCL:

* Plotting is handled through matplotlib (with cartopy or basemap), or PyNGL.
* Metadata is optional and can be disabled

wrf-python is a work in progress and should be considered "pre-alpha". Changes to the API are likely to occur between now and the first official release.

2.0 Data Structures : numpy.ndarray and xarray.DataArray¶

numpy.ndarray:¶

• Standard array object in Python for array-based computations.
• Supports all standard compiled types and also Python types.
• Has methods to reshape, perform mathematical calculations, elementwise operations, etc.
• Supports missing/fill values using the numpy.ma.MaskedArray subclass
• Does not support metadata.

xarray.DataArray (if enabled):¶

• Adds metadata capabilities to the numpy array by wrapping around it.
• HAS A numpy array, is not a numpy array subclass.
• Supports many of the numpy methods, but generally requires the numpy array to be extracted out for computations.
• Always uses NaN to represent missing/fill values, rather than the MaskedArray type familiar to numpy users.

to_np method:¶

1. If no missing/fill values, simply extracts the numpy array using the DataArray.values property.
2. If missing/fill values are found, converts the NaN values to the fill values and returns a MaskedArray.

3.0 Examples¶

Example 1: Basic Variable Extraction with getvar and Print the DataArray and Numpy Array Result¶

In [1]:

# Jupyter Notebook specific command to allow matplotlib images to be shown inline
%matplotlib inline

In [2]:

import numpy as np
import matplotlib.pyplot as plt
from Nio import open_file
from wrf import getvar, to_np

# Open the output NetCDF file with PyNIO
filename = "wrfout_d01_2010-06-13_21-00-00"
pynio_filename = filename + ".nc"
ncfile = open_file(pynio_filename)

# Alternative using netCDF4-python (for reference)
# Do 'conda install netcdf4' if the netcdf4 package is not installed
#from netCDF4 import Dataset as nc
#filename = "wrfout_d01_2010-06-13_21-00-00"
#ncfile = nc(filename)

# Extract the terrain height, which will be returned as an xarray.DataArray object (if available).  DataArray
# objects include meta data, similar to NCL's variables.
# Note:  can also use the 'ter' variable
terrainx = getvar(ncfile, "HGT", timeidx=0)

print (terrainx)

# To extract the numpy array, use the to_np function
terrain_numpy = to_np(terrainx)
print ("\nExtracted numpy array:\n")
print (terrain_numpy)

---------------------------------------------------------------------------
NIOError                                  Traceback (most recent call last)
<ipython-input-2-9e6cb63ecf8c> in <module>()
      7 filename = "wrfout_d01_2010-06-13_21-00-00"
      8 pynio_filename = filename + ".nc"
----> 9 ncfile = open_file(pynio_filename)
     10 
     11 # Alternative using netCDF4-python (for reference)

/Users/ladwig/miniconda2/lib/python2.7/site-packages/PyNIO/Nio.pyc in open_file(filename, mode, options, history, format)
    733     mask_above_value = _get_option_value(options,_Nio.option_defaults,'MaskAboveValue')
    734 
--> 735     file = _Nio.open_file(filename,mode,options,history,format)
    736 
    737     file_proxy = _proxy(file, 'str', __del__=__del__,create_variable=create_variable,create_group=create_group,close=close)

NIOError: Unable to open file

Example 2: Use getvar to Extract the 'HGT' Variable and Make a Plot Using Matplotlib with Basemap¶

In [ ]:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
from Nio import open_file
from wrf import getvar, to_np, get_basemap, latlon_coords

# Open the output NetCDF file with PyNIO
filename = "wrfout_d01_2010-06-13_21-00-00"
pynio_filename = filename + ".nc"
ncfile = open_file(pynio_filename)

# Extract the terrain height, which will be returned as an xarray.DataArray object (if available).  DataArray
# objects include meta data, similar to NCL's variables.
# Note:  can also use the 'ter' variable
terrainx = getvar(ncfile, "HGT", timeidx=0)

# Use to_np to extract the numpy array, since matplotlib does not handle xarray.DataArray natively
terrain_data = to_np(terrainx)

# Get the lat/lon 2D coordinate arrays.  Use to_np to extract the numpy array since basemap does not
# handle xarray.DataArray natively.
lons = latlon_coords(terrainx, as_np=True)
lats = latlon_coords(terrainx, as_np=True)

# Extract the basemap object from the projection information
bm = get_basemap(terrainx)

# Convert the lat/lon coordinates to projected x,y
x,y = bm(lons, lats)

# Create the figure
fig = plt.figure(figsize=(16,16))
ax = fig.add_axes([0.1,0.1,0.8,0.8])

# Draw filled contours from 100 to 3000 m, every 200 meters.
levels = np.arange(100, 3000, 200)
bm.contourf(x, y, terrain_data, levels=levels, extend="max", cmap=get_cmap("terrain"))

# Draw the coastlines and country borders.
bm.drawcoastlines()
bm.drawcountries()

# Draw the color bar
plt.colorbar(ax=ax, shrink=.7)

# Add a title
plt.title("Terrain Height (m)", {"fontsize" : 20})

plt.show()

Example 3: Use getvar to Compute Dewpoint and Make a Plot Using Matplotlib with Basemap¶

In [ ]:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
from Nio import open_file
from wrf import getvar, to_np, get_basemap, latlon_coords

# Open the output NetCDF file with PyNIO
filename = "wrfout_d01_2010-06-13_21-00-00"
pynio_filename = filename + ".nc"
ncfile = open_file(pynio_filename)

# Extract the dewpoint, which will be returned as an xarray.DataArray object (if available).  DataArray
# objects include meta data, similar to NCL's variables.
dewpointx = getvar(ncfile, "td", timeidx=0)


# Dewpoint is a 3D variable, so let's just use the lowest level
dewpointx_sfc = dewpointx[0,:,:]

# Use to_np to extract the numpy array, since matplotlib does not handle xarray.DataArray natively
dewpoint_ndarray = to_np(dewpointx_sfc)

# Get the lat/lon 2D coordinate arrays.  Use to_np to extract the numpy array since basemap does not
# handle xarray.DataArray natively.
lons = latlon_coords(dewpointx_sfc, as_np=True)
lats = latlon_coords(dewpointx_sfc, as_np=True)

# Extract the basemap object from the projection information
bm = get_basemap(dewpointx_sfc)

# Convert the lat/lon coordinates to projected x,y
x,y = bm(lons, lats)

# Create the figure
fig = plt.figure(figsize=(16,16))
ax = fig.add_axes([0.1,0.1,0.8,0.8])

# Draw filled contours from -20 C to 40 C, every 5 C.
levels = np.arange(-20, 40, 5)
bm.contourf(x, y, dewpoint_ndarray, levels=levels, extend="both", cmap=get_cmap("RdYlGn"))

# Draw the coastlines, country borders, and states.
bm.drawcoastlines()
bm.drawcountries()
bm.drawstates()

# Draw the color bar
plt.colorbar(ax=ax, shrink=.7)

# Add a title
plt.title("Dewpoint Temperature (C)", {"fontsize" : 20})

plt.show()

Example 4: Create a Vertical Cross-Section Using vertcross and Make a Plot Using Matplotlib with Basemap¶

In [ ]:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
from Nio import open_file
from wrf import getvar, vertcross, to_np

# Open the output NetCDF file with PyNIO
filename = "wrfout_d01_2010-06-13_21-00-00"
pynio_filename = filename + ".nc"
ncfile = open_file(pynio_filename)

# Extract pressure and model height
p = getvar(ncfile, "pressure", timeidx=0)
z = getvar(ncfile, "z", timeidx=0)

# Define a horizontal cross section going left to right using a pivot point in the center of the grid
# with an angle of 90.0 degrees (0 degrees is vertical).
# Pivot point is a tuple of (x, y)
pivot_point = (z.shape[-1] / 2, z.shape[-2] / 2) 
angle = 90.0

# Compute the vertical cross-section interpolation.  Also, include the lat/lon points along the cross-section.
p_vertx = vertcross(p, z, pivot_point=pivot_point, angle=angle, latlon=True)

# Extract the numpy array
p_vert_array = to_np(p_vertx)

# Create the figure
fig = plt.figure(figsize=(20,8))
ax = plt.axes([0.1,0.1,0.8,0.8])

# Define the levels [0, 50, 100, 150, ....]
levels = [0 + 50*n for n in range(20)]

# Make the contour plot
plt.contour(p_vert_array, levels=levels)

# Add the color bar
plt.colorbar(ax=ax)

# Set the x-ticks to use latitude and longitude labels.
coord_pairs = to_np(p_vertx.coords["xy_loc"])
x_ticks = np.arange(coord_pairs.shape[0])
x_labels = [pair.latlon_str() for pair in to_np(coord_pairs)]
plt.xticks(x_ticks[::100], x_labels[::100]) # Only use every 100th tick.

# Set the y-ticks to be height.
vert_vals = to_np(p_vertx.coords["vertical"])
v_ticks = np.arange(vert_vals.shape[0])
plt.yticks(v_ticks[::10], vert_vals[::10]) # Only use every 10th tick.

# Set the x-axis and  y-axis labels
plt.xlabel("Latitude, Longitude", fontsize=14)
plt.ylabel("Height (m)", fontsize=14)

# Add a title
plt.title("Vertical Cross-Section of Pressure (hPa)", {"fontsize" : 20})

plt.show()

Example 5: Interpolate Wind and Height to the 500 hPa Level Using interplevel and Make a Plot of the 500 hPa Wind Barbs, Wind Speed, and Height Contours Using Matplotlib with Basemap¶

In [ ]:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
from Nio import open_file
from wrf import getvar, interplevel, to_np, get_basemap, latlon_coords

# Open the output NetCDF file with PyNIO
filename = "wrfout_d01_2010-06-13_21-00-00"
pynio_filename = filename + ".nc"
ncfile = open_file(pynio_filename)

# Extract pressure, model height, u and v winds on mass points
p = getvar(ncfile, "pressure", timeidx=0)
z = getvar(ncfile, "z", timeidx=0, units="dm")
ua = getvar(ncfile, "ua", timeidx=0, units="kts")
va = getvar(ncfile, "va", timeidx=0, units="kts")
wspd = getvar(ncfile, "wspd_wdir", timeidx=0, units="kts")[0,...]

# Interpolate height, u, and v to to 500 hPa
ht_500 = interplevel(z, p, 500)
u_500 = interplevel(ua, p, 500)
v_500 = interplevel(va, p, 500)
wspd_500 = interplevel(wspd, p, 500)

# Get the projection
bm = get_basemap(p)

# Basemap needs numpy arrays, extract with to_np
lons = latlon_coords(ht_500, as_np=True)
lats = latlon_coords(ht_500, as_np=True)

# Convert the lat/lon coordinates to projected x,y
x,y = bm(lons, lats)

fig = plt.figure(figsize=(20,20))
ax = plt.axes([0.1,0.1,0.8,0.8])

# Draw the coastlines and country borders.
bm.drawcoastlines()
bm.drawcountries()
bm.drawstates()

# Make the 500 hPa height contours
ht_contours = bm.contour(x, y, to_np(ht_500), 10, linewidths=2.0, colors="black")

# Use contour labels for height
plt.clabel(ht_contours, inline=True, fontsize=12, fmt="%i")

# Make the wind speed filled contours
levels = np.arange(40, 120, 10)
bm.contourf(x, y, to_np(wspd_500), levels=levels, extend="max", cmap=get_cmap("rainbow"))

# Make the wind barbs.  Only use every 50th in each direction.
bm.barbs(x[::50,::50], y[::50,::50], to_np(u_500[::50, ::50]), to_np(v_500[::50, ::50]))

# Make the color bar
plt.colorbar(ax=ax, shrink=.7)

# Make the title
plt.title("500 MB Heights (dm), Wind Speed (kts), and Wind Barbs (kts)", {"fontsize" : 20})

plt.show()

In [ ]: