wrf-python/test/ipynb/WRF_Workshop_Demo.ipynb

In [ ]:

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

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, npvalues

#filename = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-02-25_18_00_00"
#filename = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2010-06-13_21:00:00"
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)

print (terrainx)

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

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, npvalues


#filename = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-02-25_18_00_00"
#filename = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2010-06-13_21:00:00"
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 npvalues to extract the numpy array, since matplotlib does not handle xarray.DataArray natively
terrain_data = npvalues(terrainx)

# Get the lat/lon 2D coordinate arrays.  Use npvalues to extract the numpy array since basemap does not
# handle xarray.DataArray natively.
lons = npvalues(terrainx.coords["XLONG"])
lats = npvalues(terrainx.coords["XLAT"])

# Extract the basemap object from the projection information
wrf_proj = terrainx.attrs["projection"]
bm = wrf_proj.basemap()

# 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()

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, npvalues


#filename = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-02-25_18_00_00"
#filename = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2010-06-13_21:00:00"
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 npvalues to extract the numpy array, since matplotlib does not handle xarray.DataArray natively
dewpoint_ndarray = npvalues(dewpointx_sfc)

# Get the lat/lon 2D coordinate arrays.  Use npvalues to extract the numpy array since basemap does not
# handle xarray.DataArray natively.
lons = npvalues(dewpointx_sfc.coords["XLONG"])
lats = npvalues(dewpointx_sfc.coords["XLAT"])

# Extract the basemap object from the projection information
wrf_proj = dewpointx_sfc.attrs["projection"]
bm = wrf_proj.basemap()

# 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()

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, npvalues

# Open the output netcdf file
#filename = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-02-25_18_00_00"
#filename = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2010-06-13_21:00:00"
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")
z = getvar(ncfile, "z")

# 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, include_latlon=True)

# Extract the numpy array
p_vert_array = npvalues(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 = npvalues(p_vertx.coords["xy_loc"])
x_ticks = np.arange(coord_pairs.shape[0])
x_labels = [pair.latlon_str() for pair in npvalues(coord_pairs)]
plt.xticks(x_ticks[::100], x_labels[::100]) # Only use every 100th tick.

# Set the y-ticks to be height.
vert_vals = npvalues(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()

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, npvalues

# Open the output netcdf file
#filename = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-02-25_18_00_00"
#filename = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2010-06-13_21:00:00"
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")
z = getvar(ncfile, "z", units="dm")
ua = getvar(ncfile, "ua", units="kts")
va = getvar(ncfile, "va", units="kts")
wspd = getvar(ncfile, "wspd_wdir", 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
wrf_proj = p.attrs["projection"]
bm = wrf_proj.basemap()

# Basemap needs numpy arrays, extract with npvalues
lons = npvalues(ht_500.coords["XLONG"])
lats = npvalues(ht_500.coords["XLAT"])

# 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, npvalues(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, npvalues(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], npvalues(u_500[::50, ::50]), npvalues(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 [ ]: