Renamed the geobounds function since it was confusing the autosummary tool with Sphinx. Updated documentation.

doc/source/basic_usage.rst

@@ -19,7 +19,7 @@ In the example below, sea level pressure is calculated and printed.
 
 .. code-block:: python
 
-    from __future__ import (absolute_import, division, print_function, unicode_literals)
+    from __future__ import print_function
 
     from netCDF4 import Dataset
     from wrf import getvar
@@ -78,7 +78,7 @@ NetCDF variables.
 
 .. code-block:: python
 
-    from __future__ import (absolute_import, division, print_function, unicode_literals)
+    from __future__ import print_function
 
     from netCDF4 import Dataset
     from wrf import getvar
@@ -141,7 +141,7 @@ The example below illustrates both.
 
 .. code-block:: python
 
-    from __future__ import (absolute_import, division, print_function, unicode_literals)
+    from __future__ import print_function
 
     from netCDF4 import Dataset
     from wrf import getvar, disable_xarray
@@ -169,7 +169,7 @@ Extracting a Numpy Array from a DataArray
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 If you need to convert an :class:`xarray.DataArray` to a :class:`numpy.ndarray`,
-wrf-python provides the :meth:`wrf.npvalues` function for this purpose. Although
+wrf-python provides the :meth:`wrf.to_np` function for this purpose. Although
 an :class:`xarray.DataArary` object already contains the 
 :attr:`xarray.DataArray.values` attribute to extract the Numpy array, there is a 
 problem when working with compiled extensions. The behavior for xarray (and pandas) 
@@ -177,18 +177,18 @@ is to convert missing/fill values to NaN, which may cause crashes when working
 with compiled extensions.  Also, some existing code may be designed to work with 
 :class:`numpy.ma.MaskedArray`, and numpy arrays with NaN may not work with it.
 
-The :meth:`wrf.npvalues` function does the following:
+The :meth:`wrf.to_np` function does the following:
 
-#. If no missing/fill values are used, :meth:`wrf.npvalues` simply returns the 
+#. If no missing/fill values are used, :meth:`wrf.to_np` simply returns the 
    :attr:`xarray.DataArray.values` attribute.
 
-#. If missing/fill values are used, then :meth:`wrf.npvalues` replaces the NaN
+#. If missing/fill values are used, then :meth:`wrf.to_np` replaces the NaN
    values with the _FillValue found in the :attr:`xarray.DataArray.attrs` 
    attribute (required) and a :class:`numpy.ma.MaskedArray` is returned.
 
 .. code-block:: python
 
-    from __future__ import (absolute_import, division, print_function, unicode_literals)
+    from __future__ import print_function
 
     from netCDF4 import Dataset
     from wrf import getvar
@@ -198,7 +198,7 @@ The :meth:`wrf.npvalues` function does the following:
     # Get the Sea Level Pressure
     cape_3d = getvar(ncfile, "cape_3d")
     
-    cape_3d_ndarray = npvalues(cape_3d)
+    cape_3d_ndarray = to_np(cape_3d)
     
     print(type(cape_3d_ndarray))
 
@@ -229,12 +229,11 @@ function.
 
 .. code-block:: python
 
-    from __future__ import (absolute_import, division, print_function, unicode_literals)
+    from __future__ import print_function
     
     from netCDF4 import Dataset
     from wrf import getvar, ALL_TIMES
     
-    
     # Creating a simple test list with three timesteps
     wrflist = [Dataset("wrfout_d01_2016-10-07_00_00_00"), 
                Dataset("wrfout_d01_2016-10-07_01_00_00"), 
@@ -305,7 +304,7 @@ for most cases.
 
 .. code-block:: python
 
-    from __future__ import (absolute_import, division, print_function, unicode_literals)
+    from __future__ import print_function
 
     from netCDF4 import Dataset
     from wrf import getvar, ALL_TIMES
@@ -367,7 +366,7 @@ numpy's automatic squeezing of the single 'Time' dimension. To maintain the
 
 .. code-block:: python
 
-    from __future__ import (absolute_import, division, print_function, unicode_literals)
+    from __future__ import print_function
 
     from netCDF4 import Dataset
     from wrf import getvar, ALL_TIMES
@@ -437,7 +436,7 @@ will be combined.
 
 .. code-block:: python
     
-    from __future__ import (absolute_import, division, print_function, unicode_literals)
+    from __future__ import print_function
 
     from netCDF4 import Dataset
     from wrf import getvar, ALL_TIMES
@@ -504,7 +503,7 @@ a specific horizontal level, usually pressure or height.
 
 .. code-block:: python
     
-    from __future__ import (absolute_import, division, print_function, unicode_literals)
+    from __future__ import print_function
     
     from netCDF4 import Dataset
     from wrf import getvar, interplevel
@@ -579,7 +578,7 @@ Example Using Start Point and End Point
 
 .. code-block:: python
 
-    from __future__ import (absolute_import, division, print_function, unicode_literals)
+    from __future__ import print_function
 
     from netCDF4 import Dataset
     from wrf import getvar, vertcross, CoordPair
@@ -644,7 +643,7 @@ Example Using Pivot Point and Angle
 
 .. code-block:: python
 
-    from __future__ import (absolute_import, division, print_function, unicode_literals)
+    from __future__ import print_function
 
     from netCDF4 import Dataset
     from wrf import getvar, vertcross, CoordPair  
@@ -710,7 +709,7 @@ Example Using Lat/Lon Coordinates
 
 .. code-block:: python
 
-    from __future__ import (absolute_import, division, print_function, unicode_literals)
+    from __future__ import print_function
 
     from netCDF4 import Dataset
     from wrf import getvar, vertcross, CoordPair  
@@ -781,7 +780,7 @@ Example Using Specified Vertical Levels
 
 .. code-block:: python
 
-    from __future__ import (absolute_import, division, print_function, unicode_literals)
+    from __future__ import print_function
 
     from netCDF4 import Dataset
     from wrf import getvar, vertcross, CoordPair  
@@ -861,7 +860,7 @@ Example Using Start Point and End Point
 
 .. code-block:: python
 
-    from __future__ import (absolute_import, division, print_function, unicode_literals)
+    from __future__ import print_function
 
     from netCDF4 import Dataset
     from wrf import getvar, interpline, CoordPair  
@@ -912,7 +911,7 @@ Example Using Pivot Point and Angle
 
 .. code-block:: python
 
-    from __future__ import (absolute_import, division, print_function, unicode_literals)
+    from __future__ import print_function
 
     from netCDF4 import Dataset
     from wrf import getvar, interpline, CoordPair  
@@ -962,7 +961,7 @@ Example Using Lat/Lon Coordinates
 
 .. code-block:: python
 
-    from __future__ import (absolute_import, division, print_function, unicode_literals)
+    from __future__ import print_function
 
     from netCDF4 import Dataset
     from wrf import getvar, interpline, CoordPair  
@@ -1024,7 +1023,7 @@ The surface levels to interpolate also need to be specified.
 
 .. code-block:: python
 
-    from __future__ import (absolute_import, division, print_function, unicode_literals)
+    from __future__ import print_function
 
     from netCDF4 import Dataset
     from wrf import getvar, interpline, CoordPair  
@@ -1097,7 +1096,7 @@ Example With Single Coordinates
 
 .. code-block:: python
 
-    from __future__ import (absolute_import, division, print_function, unicode_literals)
+    from __future__ import print_function
 
     from netCDF4 import Dataset
     from wrf import getvar, interpline, CoordPair, xy_to_ll, ll_to_xy
@@ -1137,7 +1136,7 @@ Example With Multiple Coordinates
 
 .. code-block:: python
 
-    from __future__ import (absolute_import, division, print_function, unicode_literals)
+    from __future__ import print_function
     
     from netCDF4 import Dataset
     from wrf import getvar, interpline, CoordPair, xy_to_ll, ll_to_xy

doc/source/conf.py

@@ -25,7 +25,7 @@ class Mock(MagicMock):
     def __getattr__(cls, name):
         return Mock()
     
-MOCK_MODULES = ["numpy", "numpy.ma", "xarray", 
+MOCK_MODULES = ["numpy", "numpy.ma", "xarray", "cartopy", 
                 "pandas", "matplotlib", "netCDF4", "mpl_toolkits.basemap",
                 "wrf._wrffortran"]
 sys.modules.update((mod_name, Mock()) for mod_name in MOCK_MODULES)

doc/source/new.rst

@@ -4,7 +4,17 @@ What's New
 v1.0a3
 -----------
 
-- Alpha release 3
+- Alpha release 3.
-- Added docstrings
+- Added docstrings.
-- Now uses CoordPair for cross sections so that lat/lon can be used
+- The mapping API has changed.
-- Renamed some functions and or arguments
+    - The projection attributes are no longer arrays for moving domains.
+    - Utility functions have been added for extracting geobounds.  It is now 
+      easier to get map projection objects from sliced variables.
+    - Utility functions have been added for getting cartopy, basemap, and pyngl
+      objects.
+    - Users should no longer need to use xarray attributes directly
+- Now uses CoordPair for cross sections so that lat/lon can be used instead of 
+  raw x,y grid coordinates.
+- Renamed npvalues to to_np which is more intuitive.
+- Fixed issue with generator expressions.
+- Renamed some functions and arguments.

doc/source/plot.rst

@@ -13,7 +13,7 @@ Matplotlib With Cartopy
 Cartopy is becoming the main tool for base mapping with matplotlib, but you should 
 be aware of a few shortcomings when working with WRF data.
 
-- The builtin tranformations of coordinates when calling the contouring functions
+- The builtin transformations of coordinates when calling the contouring functions
   do not work correctly with the rotated pole projection.  The 
   transform_points method needs to be called manually on the latitude and 
   longitude arrays.
@@ -34,15 +34,14 @@ Plotting a Two-dimensional Field
    
 .. code-block:: python
     
-    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
+    from wrf import (to_np, getvar, smooth2d, get_cartopy, latlon_coords,
+                    cartopy_xlim, cartopy_ylim)
     
     ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00")
     
@@ -52,15 +51,11 @@ Plotting a Two-dimensional Field
     # Smooth the sea level pressure since it tends to be noisy near the mountains
     smooth_slp = smooth2d(slp, 3)
     
-    # Get the numpy array from the XLAT and XLONG coordinates
+    # Get the latitude and longitude coordinate arrays
-    lats = npvalues(slp.coords["XLAT"])
+    lats, lons = latlon_coords(smooth_slp)
-    lons = npvalues(slp.coords["XLONG"])
-    
-    # Get the wrf.WrfProj object
-    wrf_proj = slp.attrs["projection"]
     
     # The WrfProj.cartopy() method returns a cartopy.crs projection object
-    cart_proj = wrf_proj.cartopy()
+    cart_proj = get_cartopy(smooth_slp)
     
     # Create a figure that's 10x10
     fig = plt.figure(figsize=(10,10))
@@ -75,17 +70,583 @@ Plotting a Two-dimensional Field
     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 
+    plt.contour(to_np(lons), to_np(lats), to_np(smooth_slp), 10, colors="black", 
-    # cartopy to transform the points in to the WRF projection set for the GeoAxes.
+                transform=crs.PlateCarree())
-    plt.contour(lons, lats, npvalues(smooth_slp), 10, colors="black", transform=crs.PlateCarree())
+    plt.contourf(to_np(lons), to_np(lats), to_np(smooth_slp), 10, 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_xlim(cartopy_xlim(smooth_slp))
-    ax.set_ylim(wrf_proj.cartopy_ylim())
+    ax.set_ylim(cartopy_ylim(smooth_slp))
     
     # Add the gridlines
     ax.gridlines()
+    
+    plt.title("Sea Level Pressure (hPa)")
+    
+    plt.show()
+
+Horizontal Interpolation to a Pressure Level
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. image:: _static/images/cartopy_500.png    
+   :scale: 100%
+   :align: center
+    
+.. code-block:: python
+
+    from netCDF4 import Dataset 
+    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 wrf import (getvar, interplevel, to_np, get_cartopy, 
+                    latlon_coords, cartopy_xlim, cartopy_ylim)
+    
+    ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00")
+    
+    # Extract the pressure, geopotential height, and wind variables
+    p = getvar(ncfile, "pressure")
+    z = getvar(ncfile, "z", units="dm")
+    ua = getvar(ncfile, "ua", units="kt")
+    va = getvar(ncfile, "va", units="kt")
+    wspd = getvar(ncfile, "wspd_wdir", units="kts")[0,:]
+    
+    # Interpolate geopotential height, u, and v winds 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 lat/lon coordinates
+    lats, lons = latlon_coords(ht_500)
+    
+    # Get the cartopy map projection information
+    cart_proj = get_cartopy(ht_500)
+    
+    # Create the figure
+    fig = plt.figure(figsize=(10,10))
+    ax = plt.axes([0.1,0.1,0.8,0.8], 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=0.5)
+    ax.coastlines('50m', linewidth=0.8)
+    
+    # Create the 500 hPa geopotential height contours
+    levels = np.arange(520., 580., 6.)
+    contours = plt.contour(to_np(lons), to_np(lats), to_np(ht_500), levels=levels, colors="black", 
+                transform=crs.PlateCarree())
+    plt.clabel(contours, inline=1, fontsize=10, fmt="%i")
+    
+    # Wind Speed contours
+    levels = [25, 30, 35, 40, 50, 60, 70, 80, 90, 100, 110, 120]
+    wspd_contours = plt.contourf(to_np(lons), to_np(lats), to_np(wspd_500), levels=levels,
+                                 cmap=get_cmap("rainbow"), 
+                                 transform=crs.PlateCarree())
+    plt.colorbar(wspd_contours, ax=ax, orientation="horizontal", pad=.05)
+    
+    # Add the 500 hPa wind barbs, only plotting every 125th data point.
+    plt.barbs(to_np(lons[::125,::125]), to_np(lats[::125,::125]), to_np(u_500[::125, ::125]), 
+              to_np(v_500[::125, ::125]), transform=crs.PlateCarree(), length=6)
+    
+    # Set the map bounds
+    ax.set_xlim(cartopy_xlim(ht_500))
+    ax.set_ylim(cartopy_ylim(ht_500))
+    
+    ax.gridlines()
+    
+    plt.title("500 MB Height (dm), Wind Speed (kt), Barbs (kt)")
+    
+    plt.show()
+    
+
+Panel Plots From Front Page
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+This lengthy example shows how to make the panel plots on the first page 
+of the documentation.  For a simpler example of how to make a cross section 
+plot, see :ref:`cross_example`.
+
+.. image:: _static/images/matthew.png    
+   :scale: 100%
+   :align: center
+
+.. code-block:: python
+
+    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
+    
+    from wrf import (getvar, to_np, vertcross, smooth2d, CoordPair, GeoBounds, 
+                    latlon_coords, get_cartopy, cartopy_xlim, cartopy_ylim)
+    
+    # Open the output NetCDF file
+    filename = "wrfout_d01_2016-10-07_00_00_00"
+    ncfile = DataSet(filename)
+    
+    # Get the WRF variables
+    slp = getvar(ncfile, "slp")
+    smooth_slp = smooth2d(slp, 3)
+    ctt = getvar(ncfile, "ctt")
+    height = getvar(ncfile, "z")
+    dbz = getvar(ncfile, "dbz")
+    Z = 10**(dbz/10.)
+    wspd =  getvar(ncfile, "wspd_wdir", units="kt")[0,:]
+    
+    # Set the start point and end point for the cross section
+    start_point = CoordPair(lat=26.76, lon=-80.0)
+    end_point = CoordPair(lat=26.76, lon=-77.8)
+    
+    # Compute the vertical cross-section interpolation.  Also, include the 
+    # lat/lon points along the cross-section in the metadata by setting 
+    # latlon to True.
+    z_cross = vertcross(Z, height, wrfin=ncfile, start_point=start_point, 
+                        end_point=end_point, latlon=True, meta=True)
+    wspd_cross = vertcross(wspd, height, wrfin=ncfile, start_point=start_point, 
+                           end_point=end_point, latlon=True, meta=True)
+    dbz_cross = 10.0 * np.log10(z_cross)
+    
+    # Get the latitude and longitude coordinate arrays
+    lats, lons = latlon_coords(slp)
+    
+    # Get the cartopy projection object
+    cart_proj = get_cartopy(slp)
+    
+    # Create the figure which will have 3 subplots, one on the left, and 
+    # two vertically stacked on the right.
+    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)
+    
+    # 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 ctt contours.
+    contour_levels = [960, 965, 970, 975, 980, 990]
+    c1 = ax_ctt.contour(to_np(lons), to_np(lats), to_np(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, to_np(ctt), contour_levels, 
+                                   cmap=get_cmap("Greys"), transform=crs.PlateCarree(), 
+                                   zorder=2)
+    
+    # Draw the yellow cross section line
+    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
+    hur_bounds = GeoBounds(CoordPair(lat=np.amin(to_np(lats)), lon=-85.0),
+                           CoordPair(lat=30.0, lon=-72.0))
+    ax_ctt.set_xlim(cartopy_xlim(ctt, geobounds=hur_bounds))
+    ax_ctt.set_ylim(cartopy_ylim(ctt, geobounds=hur_bounds))
+    ax_ctt.gridlines()
+    
+    # Make the contour plot for wspd
+    wspd_contours = ax_wspd.contourf(to_np(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(to_np(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 = to_np(dbz_cross.coords["xy_loc"])
+    x_ticks = np.arange(coord_pairs.shape[0])
+    x_labels = [pair.latlon_str() for pair in to_np(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 = to_np(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 titles
+    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.show()
+    
+Matplotlib with Basemap
+-----------------------
+
+Although basemap is in maintenance mode only and becoming deprecated, it is still 
+widely used by many programmers.  Cartopy is becoming the preferred package for 
+mapping, however it suffers from growing pains in some areas 
+(can't use latitude/longitude labels for many map projections).  If you 
+run in to these issues, basemap is likely to accomplish what you need, despite 
+slower performance.
+
+
+Plotting a Two-Dimensional Field
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. image:: _static/images/basemap_slp.png    
+   :scale: 100%
+   :align: center
+
+.. code-block:: python
+    
+    from netCDF4 import Dataset   
+    import matplotlib.pyplot as plt
+    from matplotlib.cm import get_cmap
+    from mpl_toolkits.basemap import Basemap
+    
+    from wrf import to_np, getvar, smooth2d, latlon_coords, get_basemap
+    
+    ncfile = Dataset("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 latitude and longitude coordinates
+    lats, lons = latlon_coords(slp)
+    
+    # Get the basemap projection object
+    bm = get_basemap(slp)
+    
+    # Create a figure that's 10x10
+    fig = plt.figure(figsize=(10,10))
+    
+    bm.drawcoastlines(linewidth=0.25)
+    bm.drawstates(linewidth=0.25)
+    bm.drawcountries(linewidth=0.25)
+    
+    # Convert the lats and lons to x and y.  Make sure you convert the lats and lons to 
+    # numpy arrays via to_np, or basemap crashes with an undefined RuntimeError.
+    x, y = bm(to_np(lons), to_np(lats))
+    
+    # Draw the contours and filled contours
+    bm.contour(x, y, to_np(smooth_slp), 10, colors="black")
+    bm.contourf(x, y, to_np(smooth_slp), 10)
+    
+    # Add a color bar
+    plt.colorbar(shrink=.47)
+    
+    plt.title("Sea Level Pressure (hPa)")
+    plt.show()
+    
+
+Horizontal Interpolation to a Pressure Level
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. image:: _static/images/basemap_500.png    
+   :scale: 100%
+   :align: center
+
+.. code-block:: python
+
+    from netCDF4 import Dataset 
+    import numpy as np
+    import matplotlib.pyplot as plt
+    from matplotlib.cm import get_cmap
+    
+    from wrf import getvar, interplevel, to_np, get_basemap, latlon_coords
+    
+    ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00")
+    
+    # Extract the pressure, geopotential height, and wind variables
+    p = getvar(ncfile, "pressure")
+    z = getvar(ncfile, "z", units="dm")
+    ua = getvar(ncfile, "ua", units="kt")
+    va = getvar(ncfile, "va", units="kt")
+    wspd = getvar(ncfile, "wspd_wdir", units="kts")[0,:]
+    
+    # Interpolate geopotential height, u, and v winds 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 lat/lon coordinates
+    lats, lons = latlon_coords(ht_500)
+    
+    # Get the basemap projection object
+    bm = get_basemap(ht_500)
+    
+    # Create the figure
+    fig = plt.figure(figsize=(8,8))
+    ax = plt.axes([0.1,0.1,0.8,0.8])
+    
+    # Get the x and y coordinates
+    x, y = bm(to_np(lons), to_np(lats))
+    
+    # Create the 500 hPa geopotential height contours
+    levels = np.arange(520., 580., 6.)
+    contours = bm.contour(x, y, to_np(ht_500), levels=levels, colors="black")
+    plt.clabel(contours, inline=1, fontsize=10, fmt="%i")
+    
+    # Wind Speed contours
+    levels = [25, 30, 35, 40, 50, 60, 70, 80, 90, 100, 110, 120]
+    wspd_contours = bm.contourf(x, y, to_np(wspd_500), levels=levels,
+                                 cmap=get_cmap("rainbow"))
+    plt.colorbar(wspd_contours, ax=ax, orientation="horizontal", pad=.05)
+    
+    bm.drawcoastlines(linewidth=0.25)
+    bm.drawstates(linewidth=0.25)
+    bm.drawcountries(linewidth=0.25)
+    
+    # Add the 500 hPa wind barbs, only plotting every 125th data point.
+    bm.barbs(x[::125,::125], y[::125,::125], to_np(u_500[::125, ::125]), 
+              to_np(v_500[::125, ::125]), length=6)
+    
+    plt.title("500 MB Height (dm), Wind Speed (kt), Barbs (kt)")
+
+    plt.show()
+    
+Panel Plots from the Front Page
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+This lengthy example shows how to make the panel plots on the first page 
+of the documentation.  For a simpler example of how to make a cross section 
+plot, see :ref:`cross_example`.
+
+.. image:: _static/images/basemap_front.png    
+   :scale: 100%
+   :align: center
+
+.. code-block:: python
+
+    import numpy as np
+    import matplotlib.pyplot as plt
+    from matplotlib.cm import get_cmap
+    from netCDF4 import Dataset
+    
+    from wrf import getvar, to_np, vertcross, smooth2d, CoordPair, get_basemap, latlon_coords
+    
+    # Open the output NetCDF file
+    filename = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-10-07_00_00_00"
+    ncfile = Dataset(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, "wspd_wdir", units="kt")[0,:]
+    
+    # Set the start point and end point for the cross section
+    start_point = CoordPair(lat=26.76, lon=-80.0)
+    end_point = CoordPair(lat=26.76, lon=-77.8)
+    
+    # Compute the vertical cross-section interpolations.  Also, include the lat/lon points along the cross-section 
+    # in the metadata by setting latlon to True.
+    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
+    lats, lons = latlon_coords(slp)
+    
+    # Create a figure that will have 3 subplots.  One on the left, 
+    # two vertically stacked on the right.
+    fig = plt.figure(figsize=(8,5))
+    ax_ctt = fig.add_subplot(1,2,1)
+    ax_wspd = fig.add_subplot(2,2,2)
+    ax_dbz = fig.add_subplot(2,2,4)
+    
+    # Get the basemap projection class
+    bm = get_basemap(slp)
+    
+    # Get the x, y values
+    x, y = bm(to_np(lons), to_np(lats))
+    
+    # Make the pressure contours.
+    contour_levels = [960, 965, 970, 975, 980, 990]
+    c1 = bm.contour(x, y, to_np(smooth_slp), levels=contour_levels, colors="white", 
+                    zorder=3, linewidths=1.0, ax=ax_ctt)
+    
+    # Create the filled cloud top temperature contours
+    contour_levels = [-80.0, -70.0, -60, -50, -40, -30, -20, -10, 0, 10]
+    ctt_contours = bm.contourf(x, y, to_np(ctt), contour_levels, cmap=get_cmap("Greys"),
+                                   zorder=2, ax=ax_ctt)
+    
+    point_x, point_y = bm([start_point.lon, end_point.lon], [start_point.lat, end_point.lat])
+    bm.plot([point_x[0], point_x[1]], [point_y[0], point_y[1]], color="yellow", 
+                marker="o", zorder=3, ax=ax_ctt)
+    
+    # Create the label bar for cloud top temperature
+    cb_ctt = fig.colorbar(ctt_contours, ax=ax_ctt, shrink=.68)
+    cb_ctt.ax.tick_params(labelsize=5)
+    
+    # Draw the oceans, land, and states.  Use the same colors as the cartopy
+    # example
+    bm.drawcoastlines(linewidth=0.25, ax=ax_ctt)
+    bm.drawstates(linewidth=0.25, ax=ax_ctt)
+    bm.drawcountries(linewidth=0.25, ax=ax_ctt)
+    bm.fillcontinents(color=np.array([ 0.9375 , 0.9375 , 0.859375]), 
+                      ax=ax_ctt, lake_color=np.array([ 0.59375 , 0.71484375, 0.8828125 ]))
+    bm.drawmapboundary(fill_color=np.array([ 0.59375 , 0.71484375, 0.8828125 ]), ax=ax_ctt)
+    
+    # Draw Parallels
+    parallels = np.arange(np.amin(lats), 30., 2.5)
+    bm.drawparallels(parallels, ax=ax_ctt)
+    
+    merids = np.arange(-85.0, -72.0, 2.5)
+    bm.drawmeridians(merids, ax=ax_ctt)
+    
+    # Get the x and y coordinate ranges for the cropped region around the 
+    # hurricane
+    x_start, y_start = bm(-85.0, np.amin(lats))
+    x_end, y_end = bm(-72.0, 30.0)
+    
+    ax_ctt.set_xlim([x_start, x_end])
+    ax_ctt.set_ylim([y_start, y_end])
+    
+    # Make the contour plot for wspd
+    wspd_contours = ax_wspd.contourf(to_np(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(to_np(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 = to_np(dbz_cross.coords["xy_loc"])
+    x_ticks = np.arange(coord_pairs.shape[0])
+    x_labels = [pair.latlon_str() for pair in to_np(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 = to_np(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.show()
+    
+    
+.. _cross_example:
+
+Vertical Cross Section
+-------------------------------
+
+Vertical cross sections require no mapping software package and can be 
+plotted using the standard matplotlib API.
+
+.. image:: _static/images/cartopy_cross.png    
+   :scale: 100%
+   :align: center
+
+.. code-block:: python
+
+    import numpy as np
+    import matplotlib.pyplot as plt
+    from matplotlib.cm import get_cmap
+    from netCDF4 import Dataset
+    
+    from wrf import to_np, getvar, CoordPair, vertcross
+    
+    # Open the output NetCDF file with PyNIO
+    filename = "wrfout_d01_2016-10-07_00_00_00"
+    ncfile = Dataset(filename)
+    
+    # Extract pressure and model height
+    z = getvar(ncfile, "z")
+    wspd =  getvar(ncfile, "uvmet_wspd_wdir", units="kt")[0,:]
+    
+    # Define the start point and end point for the cross section using 
+    # latitude and longitude coordinates.
+    start_point = CoordPair(lat=26.76, lon=-80.0)
+    end_point = CoordPair(lat=26.76, lon=-77.8)
+    
+    # Compute the vertical cross-section interpolation.  Also, include the lat/lon 
+    # points along the cross-section in the metadata by setting latlon to True.
+    wspd_cross = vertcross(wspd, z, wrfin=ncfile, start_point=start_point, 
+                           end_point=end_point, latlon=True, meta=True)
+    
+    # Create the figure
+    fig = plt.figure(figsize=(10,5))
+    ax = plt.axes([0.1,0.1,0.8,0.8])
+    
+    # Make the contour plot
+    wspd_contours = ax.contourf(to_np(wspd_cross))
+    # Add the color bar
+    plt.colorbar(wspd_contours, ax=ax)
+    
+    # Set the x-ticks to use latitude and longitude labels.
+    coord_pairs = to_np(wspd_cross.coords["xy_loc"])
+    x_ticks = np.arange(coord_pairs.shape[0])
+    x_labels = [pair.latlon_str(fmt="{:.2f}, {:.2f}") for pair in to_np(coord_pairs)]
+    ax.set_xticks(x_ticks[::20])
+    ax.set_xticklabels(x_labels[::20], rotation=45, fontsize=8) 
+    
+    # Set the y-ticks to be height.
+    vert_vals = to_np(wspd_cross.coords["vertical"])
+    v_ticks = np.arange(vert_vals.shape[0])
+    ax.set_yticks(v_ticks[::20])
+    ax.set_yticklabels(vert_vals[::20], fontsize=8) 
+    
+    # Set the x-axis and  y-axis labels
+    ax.set_xlabel("Latitude, Longitude", fontsize=12)
+    ax.set_ylabel("Height (m)", fontsize=12)
+    
+    plt.title("Vertical Cross Section of Wind Speed (kt)")
+    
+    plt.show()
+

doc/source/user_api/index.rst

@@ -60,7 +60,7 @@ the array object to a compiled extension.
    :nosignatures:
    :toctree: ./generated/
    
-   wrf.npvalues
+   wrf.to_np
    
     
 Variable Extraction Routines
@@ -81,6 +81,22 @@ WRF NetCDF file (or a sequence of NetCDF files).
     wrf.extract_times
     
     
+Plotting Helper Routines
+^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The routines below are used to assist with plotting.
+
+.. autosummary::
+   :nosignatures:
+   :toctree: ./generated/
+   
+    wrf.geo_bounds
+    wrf.get_cartopy
+    wrf.get_basemap
+    wrf.get_pyngl
+    wrf.cartopy_xlim
+    wrf.cartopy_ylim
+    
 Raw Diagnostic Routines
 ^^^^^^^^^^^^^^^^^^^^^^^^^^
 
@@ -226,6 +242,16 @@ CoordPair Methods
    wrf.CoordPair.latlon_str
    wrf.CoordPair.xy_str
    
+GeoBounds Class
+^^^^^^^^^^^^^^^^^^^^^^^
+
+The class below is used for specifying geographic boundaries. 
+
+.. autosummary::
+   :nosignatures:
+   :toctree: ./generated/
+   
+   wrf.GeoBounds
    
 Projection Classes
 ^^^^^^^^^^^^^^^^^^^^^^^^

src/wrf/api.py

@@ -24,7 +24,7 @@ from .util import (to_np, extract_global_attrs, is_standard_wrf_var,
                    from_var, combine_dims, either, get_iterable,
                    IterWrapper, is_coordvar, latlon_coordvars, is_mapping,
                    has_time_coord, is_multi_file, is_multi_time_req,
-                   get_coord_pairs, is_time_coord_var, geobounds, 
+                   get_coord_pairs, is_time_coord_var, geo_bounds, 
                    get_cartopy, get_basemap, get_pyngl, cartopy_xlim,
                    cartopy_ylim, latlon_coords)
 from .geobnds import GeoBounds, NullGeoBounds
@@ -63,7 +63,7 @@ __all__ += ["to_np", "extract_global_attrs", "is_standard_wrf_var",
             "from_var", "combine_dims", "either", "get_iterable",
             "IterWrapper", "is_coordvar", "latlon_coordvars", "is_mapping",
             "has_time_coord", "is_multi_file", "is_multi_time_req",
-            "get_coord_pairs", "is_time_coord_var", "geobounds", 
+            "get_coord_pairs", "is_time_coord_var", "geo_bounds", 
             "get_cartopy", "get_basemap", "get_pyngl", "cartopy_xlim",
             "cartopy_ylim", "latlon_coords"]
 __all__ += ["GeoBounds", "NullGeoBounds"]

src/wrf/geobnds.py

@@ -4,6 +4,19 @@ from __future__ import (absolute_import, division, print_function,
 from .coordpair import CoordPair
 
 class GeoBounds(object):
+    """A class that stores the geographic boundaries.
+    
+    Currently, only corner points are used, specified as the bottom left and 
+    top right corners.  Users can specify the corner points directly, or 
+    specify two-dimensional latitude and longitude arrays and the corner points
+    will be extracted from them.
+        
+    Attributes:
+    
+        bottom_left (:class:`wrf.CoordPair`): The bottom left coordinate.
+        top_right (:class:`wrf.CoordPair`): The top right coordinate.
+    
+    """
     def __init__(self, bottom_left=None, top_right=None, lats=None, lons=None):
         """ Initialize a :class:`wrf.GeoBounds` object.
         
@@ -29,11 +42,27 @@ class GeoBounds(object):
         if bottom_left is not None and top_right is not None:
             self.bottom_left = bottom_left
             self.top_right = top_right
+            
+            # Make sure the users set lat/lon coordinates
+            if self.bottom_left.lat is None:
+                raise ValueError("'bottom_left' parameter does not contain a "
+                                 "'lat' attribute")
+            if self.bottom_left.lon is None:
+                raise ValueError("'bottom_left' parameter does not contain a"
+                                 "'lon' attribute")
+            if self.top_right.lat is None:
+                raise ValueError("'top_right' parameter does not contain a"
+                                 "'lat' attribute")
+            if self.top_right.lon is None:
+                raise ValueError("'top_right' parameter does not contain a"
+                                 "'lon' attribute")
         elif lats is not None and lons is not None:
             self.bottom_left = CoordPair(lat=lats[0,0], lon=lons[0,0])
             self.top_right = CoordPair(lat=lats[-1,-1], lon=lons[-1,-1])
         else:
-            raise ValueError("invalid corner point arguments")
+            raise ValueError("must specify either 'bottom_top' and "
+                             "'top_right' parameters "
+                             "or 'lats' and 'lons' parameters")
         
     def __repr__(self):
         argstr = "{}, {}".format(repr(self.bottom_left), 
@@ -43,7 +72,15 @@ class GeoBounds(object):
 
  
 class NullGeoBounds(GeoBounds):
+    """An emtpy :class:`wrf.GeoBounds` subclass.
+    
+    This is used for initializing arrays of :class:`wrf.GeoBounds`, in 
+    particular when working with moving domains and variables combined with the 
+    'join' method.
+    
+    """
     def __init__(self):
+        """ Initialize a :class:`wrf.NullGeoBounds` object."""
         pass
     
     def __repr__(self):

src/wrf/util.py

@@ -2892,7 +2892,7 @@ def get_id(obj):
     return {key : get_id(val) for key,val in viewitems(obj)}
 
 
-def geobounds(var=None, wrfin=None, varname=None, timeidx=0, method="cat",
+def geo_bounds(var=None, wrfin=None, varname=None, timeidx=0, method="cat",
               squeeze=True, cache=None):
     """Return the geographic boundaries for the variable or file(s).
     
@@ -2994,7 +2994,7 @@ def geobounds(var=None, wrfin=None, varname=None, timeidx=0, method="cat",
             if xarray_enabled():
                 var = extract_vars(wrfin, timeidx, varname, method, squeeze, 
                                    cache, meta=True, _key=_key)[varname]
-                return geobounds(var)
+                return geo_bounds(var)
             else:
                 lat_coord, lon_coord, _ = _get_coord_names(wrfin, varname)
         else:
@@ -3085,10 +3085,10 @@ def _get_wrf_proj_geobnds(var, wrfin, varname, timeidx, method, squeeze,
     """
     # Using a variable
     if var is not None:
-        geobnds = geobounds(var)
+        geobnds = geo_bounds(var)
         wrf_proj = var.attrs["projection"]
     else:
-        geobnds = geobounds(wrfin=wrfin, varname=varname, timeidx=timeidx,
+        geobnds = geo_bounds(wrfin=wrfin, varname=varname, timeidx=timeidx,
                             method=method, cache=cache)
         proj_params = get_proj_params(wrfin)
         wrf_proj = getproj(**proj_params)