nilyin
/
wrf-python
forked from 3rdparty/wrf-python
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
A collection of diagnostic and interpolation routines for use with output from the Weather Research and Forecasting (WRF-ARW) Model.
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
84 Commits
5 Branches
47 Tags
40 MiB
 
 
 
 
 
 
Tree: 29a9f70605
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from '29a9f70605'
${ noResults }
wrf-python/doc/source/plot.rst

652 lines
24 KiB
Raw Blame History

Plotting Examples
=================
The examples below show how wrf-python can be used to make plots with
matplotlib (with basemap and cartopy) and PyNGL. None of these examples
make use of xarray's builtin plotting functions, since additional work is most
likely needed to extend xarray in order to work correctly. This is planned
for a future release.
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 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.
- The rotated pole projection requires the x and y limits to be set manually
using set_xlim and set_ylim.
- You can't place latitude and longitude labels on the axes when using
any projection other than Mercator or LatLon.
Plotting a Two-dimensional Field
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. image:: _static/images/cartopy_slp.png
:scale: 100%
:align: center
.. code-block:: python
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 (to_np, getvar, smooth2d, get_cartopy, latlon_coords,
cartopy_xlim, cartopy_ylim)
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 noisy near the mountains
smooth_slp = smooth2d(slp, 3)
# Get the latitude and longitude coordinate arrays
lats, lons = latlon_coords(smooth_slp)
# The WrfProj.cartopy() method returns a cartopy.crs projection object
cart_proj = get_cartopy(smooth_slp)
# 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.
plt.contour(to_np(lons), to_np(lats), to_np(smooth_slp), 10, colors="black",
transform=crs.PlateCarree())
plt.contourf(to_np(lons), to_np(lats), to_np(smooth_slp), 10, transform=crs.PlateCarree())
# Add a color bar
plt.colorbar(ax=ax, shrink=.47)
# Set the map limits
ax.set_xlim(cartopy_xlim(smooth_slp))
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()