wrf-python/test/ipynb/Doc_Examples.ipynb

In [ ]:

%matplotlib inline

In [ ]:

from __future__ import (absolute_import, division, print_function, unicode_literals)

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
import cartopy.crs as crs
import cartopy.feature as cfeature
from netCDF4 import Dataset as nc

from wrf import npvalues, getvar, smooth2d

# Open the output NetCDF with netcdf-python
filename = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-10-07_00_00_00"
ncfile = nc(filename)

# Get sea level pressure and cloud top temperature
slp = getvar(ncfile, "slp")
ctt = getvar(ncfile, "ctt")

# Smooth the SLP
smooth_slp = smooth2d(slp, 3)

# Extract the latitude and longitude coordinate arrays as regular numpy array instead of xarray.DataArray
lons = npvalues(slp.coords["XLONG"])
lats = npvalues(slp.coords["XLAT"])

# Get the cartopy projection class
wrf_proj = slp.attrs["projection"]
cart_proj = wrf_proj.cartopy()

# Create the figure
fig = plt.figure(figsize=(4,4))
ax = plt.axes([0.1,0.1,0.8,0.8], projection=cart_proj)

# Download and create the states, land, and oceans using cartopy features
states = cfeature.NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',
                             name='admin_1_states_provinces_shp')
land = cfeature.NaturalEarthFeature(category='physical', name='land', scale='50m', 
                                    facecolor=cfeature.COLORS['land'])
ocean = cfeature.NaturalEarthFeature(category='physical', name='ocean', scale='50m', 
                                     facecolor=cfeature.COLORS['water'])

# Make the pressure contours.
contour_levels = [960, 965, 970, 975, 980, 990]
c1 = plt.contour(lons, lats, npvalues(smooth_slp), levels=contour_levels, colors="white", 
            transform=crs.PlateCarree(), zorder=3, linewidths=1.0)

# Add pressure contour labels
#plt.clabel(c1, contour_levels, inline=True, fmt='%.0f', fontsize=7)

# Create the filled cloud top temperature contours
contour_levels = [-80, -70, -60, -50, -40, -30, -20, -10, 0, 10]
plt.contourf(lons, lats, npvalues(ctt), contour_levels, cmap=get_cmap("Greys"),
             transform=crs.PlateCarree(), zorder=2)

plt.plot([-80,-77.8], [26.75,26.75], color="yellow", marker="o", transform=crs.PlateCarree(), zorder=3)

# Create the label bar for cloud top temperature
#cb2 = plt.colorbar(ax=ax, fraction=0.046, pad=0.04)
cb2 = plt.colorbar(ax=ax)

# Draw the oceans, land, and states
ax.add_feature(ocean)
ax.add_feature(land)
ax.add_feature(states, linewidth=.5, edgecolor="black")

# Crop the domain to the region around the hurricane
ax.set_extent([-85., -75.0, np.amin(lats), 30.0], crs=crs.PlateCarree())
ax.gridlines(crs=crs.PlateCarree(), draw_labels=False)

# Add the title and show the image
#plt.title("Hurricane Matthew Cloud Top Temperature (degC) ")
plt.savefig("/Users/ladwig/Documents/workspace/wrf_python/doc/source/_static/images/matthew.png",
           transparent=True, bbox_inches="tight")
plt.show()

In [ ]:

%matplotlib inline

from __future__ import (absolute_import, division, print_function, unicode_literals)

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
import cartopy.crs as crs
from cartopy.feature import NaturalEarthFeature
from Nio import open_file

from wrf import npvalues, getvar, smooth2d, ll_to_xy, CoordPair, vertcross, getproj, get_proj_params, to_xy_coords

# Open the output NetCDF file with PyNIO
filename = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-10-07_00_00_00"
pynio_filename = filename + ".nc"
ncfile = open_file(pynio_filename)

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

wspd =  getvar(ncfile, "uvmet_wspd_wdir", units="kt")[0,:]
Z = 10**(dbz/10.)

start_point = CoordPair(lat=26.75, lon=-80.0)
end_point = CoordPair(lat=26.75, lon=-77.8)

# Compute the vertical cross-section interpolation.  Also, include the lat/lon points along the cross-section.
z_cross = vertcross(Z, z, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True, meta=True)
wspd_cross = vertcross(wspd, z, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True, meta=True)
dbz_cross = 10.0 * np.log10(z_cross)

# Create the figure
fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(4,4))
#ax = plt.axes([0.1,0.1,0.8,0.8])

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

# Make the contour plot
a = axes[0].contourf(npvalues(wspd_cross))
# Add the color bar
fig.colorbar(a, ax=axes[0])

b = axes[1].contourf(npvalues(dbz_cross), levels=levels)
fig.colorbar(b, ax=axes[1])

# Set the x-ticks to use latitude and longitude labels.
coord_pairs = npvalues(dbz_cross.coords["xy_loc"])
x_ticks = np.arange(coord_pairs.shape[0])
x_labels = [pair.latlon_str() for pair in npvalues(coord_pairs)]
axes[0].set_xticks(x_ticks[::20])
axes[0].set_xticklabels([], rotation=45)
axes[1].set_xticks(x_ticks[::20])
axes[1].set_xticklabels(x_labels[::20], rotation=45, fontsize=6) 


# Set the y-ticks to be height.
vert_vals = npvalues(dbz_cross.coords["vertical"])
v_ticks = np.arange(vert_vals.shape[0])
axes[0].set_yticks(v_ticks[::20])
axes[0].set_yticklabels(vert_vals[::20], fontsize=6) 
axes[1].set_yticks(v_ticks[::20])
axes[1].set_yticklabels(vert_vals[::20], fontsize=6) 

# Set the x-axis and  y-axis labels
axes[1].set_xlabel("Latitude, Longitude", fontsize=7)
axes[0].set_ylabel("Height (m)", fontsize=7)
axes[1].set_ylabel("Height (m)", fontsize=7)

# Add a title
axes[0].set_title("Cross-Section of Wind Speed (kt)", {"fontsize" : 10})
axes[1].set_title("Cross-Section of Reflectivity (dBZ)", {"fontsize" : 10})

plt.savefig("/Users/ladwig/Documents/workspace/wrf_python/doc/source/_static/images/matthew_cross.png",
           transparent=True, bbox_inches="tight")
plt.show()

In [ ]:

%matplotlib inline

from __future__ import (absolute_import, division, print_function, unicode_literals)

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
import cartopy.crs as crs
import cartopy.feature as cfeature
from Nio import open_file

from wrf import getvar, npvalues, vertcross, smooth2d, CoordPair

# Open the output NetCDF file with PyNIO
filename = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-10-07_00_00_00"
pynio_filename = filename + ".nc"
ncfile = open_file(pynio_filename)

# Get the WRF variables
slp = getvar(ncfile, "slp")
smooth_slp = smooth2d(slp, 3)
ctt = getvar(ncfile, "ctt")
z = getvar(ncfile, "z", timeidx=0)
dbz = getvar(ncfile, "dbz", timeidx=0)
Z = 10**(dbz/10.)
wspd =  getvar(ncfile, "uvmet_wspd_wdir", units="kt")[0,:]

# Set the start point and end point for the cross section
start_point = CoordPair(lat=26.75, lon=-80.0)
end_point = CoordPair(lat=26.75, lon=-77.8)

# Compute the vertical cross-section interpolation.  Also, include the lat/lon points along the cross-section.
z_cross = vertcross(Z, z, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True, meta=True)
wspd_cross = vertcross(wspd, z, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True, meta=True)
dbz_cross = 10.0 * np.log10(z_cross)

# Extract the latitude and longitude coordinate arrays as regular numpy array instead of xarray.DataArray
lons = npvalues(slp.coords["XLONG"])
lats = npvalues(slp.coords["XLAT"])

# Get the cartopy projection class
wrf_proj = slp.attrs["projection"]
cart_proj = wrf_proj.cartopy()

# Create the figure which will have 3 subplots
fig = plt.figure(figsize=(7,5))
ax_ctt = fig.add_subplot(1,2,1,projection=cart_proj)
ax_wspd = fig.add_subplot(2,2,2)
ax_dbz = fig.add_subplot(2,2,4)

## Plot the cloud top temperature

# Download and create the states, land, and oceans using cartopy features
states = cfeature.NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',
                             name='admin_1_states_provinces_shp')
land = cfeature.NaturalEarthFeature(category='physical', name='land', scale='50m', 
                                    facecolor=cfeature.COLORS['land'])
ocean = cfeature.NaturalEarthFeature(category='physical', name='ocean', scale='50m', 
                                     facecolor=cfeature.COLORS['water'])

# Make the pressure contours.
contour_levels = [960, 965, 970, 975, 980, 990]
c1 = ax_ctt.contour(lons, lats, npvalues(smooth_slp), levels=contour_levels, colors="white", 
            transform=crs.PlateCarree(), zorder=3, linewidths=1.0)

# Create the filled cloud top temperature contours
contour_levels = [-80.0, -70.0, -60, -50, -40, -30, -20, -10, 0, 10]
ctt_contours = ax_ctt.contourf(lons, lats, npvalues(ctt), contour_levels, cmap=get_cmap("Greys"),
             transform=crs.PlateCarree(), zorder=2)

ax_ctt.plot([start_point.lon, end_point.lon], [start_point.lat, end_point.lat], color="yellow", 
            marker="o", transform=crs.PlateCarree(), zorder=3)

# Create the label bar for cloud top temperature
cb_ctt = fig.colorbar(ctt_contours, ax=ax_ctt, shrink=.5)
cb_ctt.ax.tick_params(labelsize=5)

# Draw the oceans, land, and states
ax_ctt.add_feature(ocean)
ax_ctt.add_feature(land)
ax_ctt.add_feature(states, linewidth=.5, edgecolor="black")

# Crop the domain to the region around the hurricane
ax_ctt.set_extent([-85., -75.0, np.amin(lats), 30.0], crs=crs.PlateCarree())
ax_ctt.gridlines(crs=crs.PlateCarree(), draw_labels=False)

## Plot the cross sections

# Make the contour plot for wspd
wspd_contours = ax_wspd.contourf(npvalues(wspd_cross))
# Add the color bar
cb_wspd = fig.colorbar(wspd_contours, ax=ax_wspd)
cb_wspd.ax.tick_params(labelsize=5)

# Make the contour plot for dbz
levels = [5 + 5*n for n in range(15)]
dbz_contours = ax_dbz.contourf(npvalues(dbz_cross), levels=levels)
cb_dbz = fig.colorbar(dbz_contours, ax=ax_dbz)
cb_dbz.ax.tick_params(labelsize=5)

# Set the x-ticks to use latitude and longitude labels.
coord_pairs = npvalues(dbz_cross.coords["xy_loc"])
x_ticks = np.arange(coord_pairs.shape[0])
x_labels = [pair.latlon_str() for pair in npvalues(coord_pairs)]
ax_wspd.set_xticks(x_ticks[::20])
ax_wspd.set_xticklabels([], rotation=45)
ax_dbz.set_xticks(x_ticks[::20])
ax_dbz.set_xticklabels(x_labels[::20], rotation=45, fontsize=4) 


# Set the y-ticks to be height.
vert_vals = npvalues(dbz_cross.coords["vertical"])
v_ticks = np.arange(vert_vals.shape[0])
ax_wspd.set_yticks(v_ticks[::20])
ax_wspd.set_yticklabels(vert_vals[::20], fontsize=4) 
ax_dbz.set_yticks(v_ticks[::20])
ax_dbz.set_yticklabels(vert_vals[::20], fontsize=4) 

# Set the x-axis and  y-axis labels
ax_dbz.set_xlabel("Latitude, Longitude", fontsize=5)
ax_wspd.set_ylabel("Height (m)", fontsize=5)
ax_dbz.set_ylabel("Height (m)", fontsize=5)

# Add a title
ax_ctt.set_title("Cloud Top Temperature (degC)", {"fontsize" : 7})
ax_wspd.set_title("Cross-Section of Wind Speed (kt)", {"fontsize" : 7})
ax_dbz.set_title("Cross-Section of Reflectivity (dBZ)", {"fontsize" : 7})

plt.savefig("/Users/ladwig/Documents/workspace/wrf_python/doc/source/_static/images/matthew.png",
           transparent=True, bbox_inches="tight")
plt.show()

In [ ]:

# SLP
from __future__ import (absolute_import, division, print_function, unicode_literals)
    
from netCDF4 import Dataset   
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
import cartopy.crs as crs
from cartopy.feature import NaturalEarthFeature

from wrf import npvalues, getvar, smooth2d

ncfile = Dataset("/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-10-07_00_00_00")

# Get the sea level pressure
slp = getvar(ncfile, "slp")

# Smooth the sea level pressure since it tends to be noisey near the mountains
smooth_slp = smooth2d(slp, 3)

# Get the numpy array from the XLAT and XLONG coordinates
lats = npvalues(slp.coords["XLAT"])
lons = npvalues(slp.coords["XLONG"])

# Get the wrf.WrfProj object
wrf_proj = slp.attrs["projection"]

# The cartopy() method returns a cartopy.crs projection object
cart_proj = wrf_proj.cartopy()

# Create a figure that's 10x10
fig = plt.figure(figsize=(10,10))
# Get the GeoAxes set to the projection used by WRF
ax = plt.axes(projection=cart_proj)

# Download and add the states and coastlines
states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',
                             name='admin_1_states_provinces_shp')
ax.add_feature(states, linewidth=.5)
ax.coastlines('50m', linewidth=0.8)

# Make the contour outlines and filled contours for the smoothed sea level pressure.
# The transform keyword indicates that the lats and lons arrays are lat/lon coordinates and tells 
# cartopy to transform the points in to grid space.
plt.contour(lons, lats, npvalues(smooth_slp), 10, colors="black", transform=crs.PlateCarree())
plt.contourf(lons, lats, npvalues(smooth_slp), 10, transform=crs.PlateCarree())

# Add a color bar
plt.colorbar(ax=ax, shrink=.47)

# Set the map limits
ax.set_xlim(wrf_proj.cartopy_xlim())
ax.set_ylim(wrf_proj.cartopy_ylim())

# Add the gridlines
ax.gridlines()