@ -40,9 +40,9 @@ Plotting a Two-dimensional Field
import cartopy.crs as crs import cartopy.crs as crs
from cartopy.feature import NaturalEarthFeature from cartopy.feature import NaturalEarthFeature
from wrf import (to_np, getvar, smooth2d, get_cartopy, latlon_coords, from wrf import to_np, getvar, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim
cartopy_xlim, cartopy_ylim)
# Open the NetCDF file
ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00") ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00")
# Get the sea level pressure # Get the sea level pressure
@ -51,16 +51,15 @@ Plotting a Two-dimensional Field
# Smooth the sea level pressure since it tends to be noisy near the mountains # Smooth the sea level pressure since it tends to be noisy near the mountains
smooth_slp = smooth2d(slp, 3) smooth_slp = smooth2d(slp, 3)
# Get the latitude and longitude coordinate arrays # Get the latitude and longitude points
lats, lons = latlon_coords(smooth_slp) lats, lons = latlon_coords(slp)
# The WrfProj.cartopy() method returns a cartopy.crs projection object
cart_proj = get_cartopy(smooth_slp)
# Create a figure that's 10x10 # Get the cartopy mapping object
fig = plt.figure(figsize=(10,10) ) cart_proj = get_cartopy(slp)
# Get the GeoAxes set to the projection used by WRF # Create a figure
fig = plt.figure(figsize=(12,9))
# Set the GeoAxes to the projection used by WRF
ax = plt.axes(projection=cart_proj) ax = plt.axes(projection=cart_proj)
# Download and add the states and coastlines # Download and add the states and coastlines
@ -70,19 +69,18 @@ Plotting a Two-dimensional Field
ax.coastlines('50m', linewidth=0.8) ax.coastlines('50m', linewidth=0.8)
# Make the contour outlines and filled contours for the smoothed sea level pressure. # 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", plt.contour(to_np(lons), to_np(lats), to_np(smooth_slp), 10, colors="black", transform=crs.PlateCarree())
transform=crs.PlateCarree()) plt.contourf(to_np(lons), to_np(lats), to_np(smooth_slp), 10, transform=crs.PlateCarree(), cmap=get_cmap("jet"))
plt.contourf(to_np(lons), to_np(lats), to_np(smooth_slp), 10, transform=crs.PlateCarree())
# Add a color bar # Add a color bar
plt.colorbar(ax=ax, shrink=.47 ) plt.colorbar(ax=ax, shrink=.62 )
# Set the map limits # Set the map limits. Not really necessary, but used for demonstration.
ax.set_xlim(cartopy_xlim(smooth_slp)) ax.set_xlim(cartopy_xlim(smooth_slp))
ax.set_ylim(cartopy_ylim(smooth_slp)) ax.set_ylim(cartopy_ylim(smooth_slp))
# Add the gridlines # Add the gridlines
ax.gridlines() ax.gridlines(color="black", linestyle="dotted" )
plt.title("Sea Level Pressure (hPa)") plt.title("Sea Level Pressure (hPa)")
@ -104,9 +102,9 @@ Horizontal Interpolation to a Pressure Level
import cartopy.crs as crs import cartopy.crs as crs
from cartopy.feature import NaturalEarthFeature from cartopy.feature import NaturalEarthFeature
from wrf import (getvar, interplevel, to_np, get_cartopy, from wrf import getvar, interplevel, to_np, latlon_coords, get_cartopy, cartopy_xlim, cartopy_ylim
latlon_coords, cartopy_xlim, cartopy_ylim)
# Open the NetCDF file
ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00") ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00")
# Extract the pressure, geopotential height, and wind variables # Extract the pressure, geopotential height, and wind variables
@ -125,12 +123,12 @@ Horizontal Interpolation to a Pressure Level
# Get the lat/lon coordinates # Get the lat/lon coordinates
lats, lons = latlon_coords(ht_500) lats, lons = latlon_coords(ht_500)
# Get the cartopy map projection information # Get the map projection information
cart_proj = get_cartopy(ht_500) cart_proj = get_cartopy(ht_500)
# Create the figure # Create the figure
fig = plt.figure(figsize=(10,10 )) fig = plt.figure(figsize=(12,9 ))
ax = plt.axes([0.1,0.1,0.8,0.8], projection=cart_proj) ax = plt.axes(projection=cart_proj)
# Download and add the states and coastlines # Download and add the states and coastlines
states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none', states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',
@ -138,13 +136,13 @@ Horizontal Interpolation to a Pressure Level
ax.add_feature(states, linewidth=0.5) ax.add_feature(states, linewidth=0.5)
ax.coastlines('50m', linewidth=0.8) ax.coastlines('50m', linewidth=0.8)
# Create the 500 hPa geopotential height contours # Add the 500 hPa geopotential height contours
levels = np.arange(520., 580., 6.) levels = np.arange(520., 580., 6.)
contours = plt.contour(to_np(lons), to_np(lats), to_np(ht_500), levels=levels, colors="black", contours = plt.contour(to_np(lons), to_np(lats), to_np(ht_500), levels=levels, colors="black",
transform=crs.PlateCarree()) transform=crs.PlateCarree())
plt.clabel(contours, inline=1, fontsize=10, fmt="%i") plt.clabel(contours, inline=1, fontsize=10, fmt="%i")
# Wind S peed contours # Add the wind s peed contours
levels = [25, 30, 35, 40, 50, 60, 70, 80, 90, 100, 110, 120] 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, wspd_contours = plt.contourf(to_np(lons), to_np(lats), to_np(wspd_500), levels=levels,
cmap=get_cmap("rainbow"), cmap=get_cmap("rainbow"),
@ -186,19 +184,18 @@ plot, see :ref:`cross_example`.
import cartopy.feature as cfeature import cartopy.feature as cfeature
from netCDF4 import Dataset from netCDF4 import Dataset
from wrf import (getvar, to_np, vertcross, smooth2d, CoordPair, GeoBounds, from wrf import (getvar, to_np, vertcross, smooth2d, CoordPair, GeoBounds, get_cartopy,
latlon_coords, get_cartopy , cartopy_xlim, cartopy_ylim) latlon_coords, cartopy_xlim, cartopy_ylim)
# Open the output NetCDF file # Open the NetCDF file
filename = "wrfout_d01_2016-10-07_00_00_00" ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00")
ncfile = DataSet(filename)
# Get the WRF variables # Get the WRF variables
slp = getvar(ncfile, "slp") slp = getvar(ncfile, "slp")
smooth_slp = smooth2d(slp, 3) smooth_slp = smooth2d(slp, 3)
ctt = getvar(ncfile, "ctt") ctt = getvar(ncfile, "ctt")
height = getvar(ncfile, "z" ) z = getvar(ncfile, "z", timeidx=0 )
dbz = getvar(ncfile, "dbz") dbz = getvar(ncfile, "dbz", timeidx=0 )
Z = 10**(dbz/10.) Z = 10**(dbz/10.)
wspd = getvar(ncfile, "wspd_wdir", units="kt")[0,:] wspd = getvar(ncfile, "wspd_wdir", units="kt")[0,:]
@ -206,24 +203,20 @@ plot, see :ref:`cross_example`.
start_point = CoordPair(lat=26.76, lon=-80.0) start_point = CoordPair(lat=26.76, lon=-80.0)
end_point = CoordPair(lat=26.76, lon=-77.8) end_point = CoordPair(lat=26.76, lon=-77.8)
# Compute the vertical cross-section interpolation. Also, include the # Compute the vertical cross-section interpolation. Also, include the lat/lon points along the cross-section
# lat/lon points along the cross-section in the metadata by setting # in the metadata by setting latlon to True.
# latlon to True. z_cross = vertcross(Z, z, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True, meta=True)
z_cross = vertcross(Z, height, wrfin=ncfile, start_point=start_point, wspd_cross = vertcross(wspd, z, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True, meta=True)
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) dbz_cross = 10.0 * np.log10(z_cross)
# Get the latitude and longitude coordinate array s # Get the lat/lon points
lats, lons = latlon_coords(slp) lats, lons = latlon_coords(slp)
# Get the cartopy projection object # Get the cartopy projection object
cart_proj = get_cartopy(slp) cart_proj = get_cartopy(slp)
# Create the figure which will have 3 subplots, one on the left, and # Create a figure that will have 3 subplots
# two vertically stacked on the right. fig = plt.figure(figsize=(10,7))
fig = plt.figure(figsize=(7,5))
ax_ctt = fig.add_subplot(1,2,1,projection=cart_proj) ax_ctt = fig.add_subplot(1,2,1,projection=cart_proj)
ax_wspd = fig.add_subplot(2,2,2) ax_wspd = fig.add_subplot(2,2,2)
ax_dbz = fig.add_subplot(2,2,4) ax_dbz = fig.add_subplot(2,2,4)
@ -236,23 +229,21 @@ plot, see :ref:`cross_example`.
ocean = cfeature.NaturalEarthFeature(category='physical', name='ocean', scale='50m', ocean = cfeature.NaturalEarthFeature(category='physical', name='ocean', scale='50m',
facecolor=cfeature.COLORS['water']) facecolor=cfeature.COLORS['water'])
# Make the ctt contours. # Make the pressure contours
contour_levels = [960, 965, 970, 975, 980, 990] 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, c1 = ax_ctt.contour(lons, lats, to_np(smooth_slp), levels=contour_levels, colors="white" ,
colors="white", transform=crs.PlateCarree(), zorder=3, linewidths=1.0) transform=crs.PlateCarree(), zorder=3, linewidths=1.0)
# Create the filled cloud top temperature contours # Create the filled cloud top temperature contours
contour_levels = [-80.0, -70.0, -60, -50, -40, -30, -20, -10, 0, 10] 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, ctt_contours = ax_ctt.contourf(to_np(lons), to_np(lats), to_np(ctt), contour_levels, cmap=get_cmap("Greys"),
cmap=get_cmap("Greys"), transform=crs.PlateCarree(), transform=crs.PlateCarree(), zorder=2)
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",
ax_ctt.plot([start_point.lon, end_point.lon], [start_point.lat, end_point.lat], marker="o", transform=crs.PlateCarree(), zorder=3)
color="yellow", marker="o", transform=crs.PlateCarree(), zorder=3)
# Create the label bar for cloud top temperature # Create the color bar for cloud top temperature
cb_ctt = fig.colorbar(ctt_contours, ax=ax_ctt, shrink=.5 ) cb_ctt = fig.colorbar(ctt_contours, ax=ax_ctt, shrink=.60 )
cb_ctt.ax.tick_params(labelsize=5) cb_ctt.ax.tick_params(labelsize=5)
# Draw the oceans, land, and states # Draw the oceans, land, and states
@ -265,21 +256,21 @@ plot, see :ref:`cross_example`.
CoordPair(lat=30.0, lon=-72.0)) CoordPair(lat=30.0, lon=-72.0))
ax_ctt.set_xlim(cartopy_xlim(ctt, geobounds=hur_bounds)) ax_ctt.set_xlim(cartopy_xlim(ctt, geobounds=hur_bounds))
ax_ctt.set_ylim(cartopy_ylim(ctt, geobounds=hur_bounds)) ax_ctt.set_ylim(cartopy_ylim(ctt, geobounds=hur_bounds))
ax_ctt.gridlines() ax_ctt.gridlines(color="white", linestyle="dotted" )
# Make the contour plot for wspd # Make the contour plot for wind spee d
wspd_contours = ax_wspd.contourf(to_np(wspd_cross)) wspd_contours = ax_wspd.contourf(to_np(wspd_cross), cmap=get_cmap("jet") )
# Add the color bar # Add the color bar
cb_wspd = fig.colorbar(wspd_contours, ax=ax_wspd) cb_wspd = fig.colorbar(wspd_contours, ax=ax_wspd)
cb_wspd.ax.tick_params(labelsize=5) cb_wspd.ax.tick_params(labelsize=5)
# Make the contour plot for dbz # Make the contour plot for dbz
levels = [5 + 5*n for n in range(15)] levels = [5 + 5*n for n in range(15)]
dbz_contours = ax_dbz.contourf(to_np(dbz_cross), levels=levels) dbz_contours = ax_dbz.contourf(to_np(dbz_cross), levels=levels, cmap=get_cmap("jet") )
cb_dbz = fig.colorbar(dbz_contours, ax=ax_dbz) cb_dbz = fig.colorbar(dbz_contours, ax=ax_dbz)
cb_dbz.ax.tick_params(labelsize=5) cb_dbz.ax.tick_params(labelsize=5)
# Set the x-ticks to use latitude and longitude labels. # Set the x-ticks to use latitude and longitude labels
coord_pairs = to_np(dbz_cross.coords["xy_loc"]) coord_pairs = to_np(dbz_cross.coords["xy_loc"])
x_ticks = np.arange(coord_pairs.shape[0]) x_ticks = np.arange(coord_pairs.shape[0])
x_labels = [pair.latlon_str() for pair in to_np(coord_pairs)] x_labels = [pair.latlon_str() for pair in to_np(coord_pairs)]
@ -288,7 +279,7 @@ plot, see :ref:`cross_example`.
ax_dbz.set_xticks(x_ticks[::20]) ax_dbz.set_xticks(x_ticks[::20])
ax_dbz.set_xticklabels(x_labels[::20], rotation=45, fontsize=4) ax_dbz.set_xticklabels(x_labels[::20], rotation=45, fontsize=4)
# Set the y-ticks to be height. # Set the y-ticks to be height
vert_vals = to_np(dbz_cross.coords["vertical"]) vert_vals = to_np(dbz_cross.coords["vertical"])
v_ticks = np.arange(vert_vals.shape[0]) v_ticks = np.arange(vert_vals.shape[0])
ax_wspd.set_yticks(v_ticks[::20]) ax_wspd.set_yticks(v_ticks[::20])
@ -301,7 +292,7 @@ plot, see :ref:`cross_example`.
ax_wspd.set_ylabel("Height (m)", fontsize=5) ax_wspd.set_ylabel("Height (m)", fontsize=5)
ax_dbz.set_ylabel("Height (m)", fontsize=5) ax_dbz.set_ylabel("Height (m)", fontsize=5)
# Add titles # Add a title
ax_ctt.set_title("Cloud Top Temperature (degC)", {"fontsize" : 7}) ax_ctt.set_title("Cloud Top Temperature (degC)", {"fontsize" : 7})
ax_wspd.set_title("Cross-Section of Wind Speed (kt)", {"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}) ax_dbz.set_title("Cross-Section of Reflectivity (dBZ)", {"fontsize" : 7})
@ -333,25 +324,27 @@ Plotting a Two-Dimensional Field
from matplotlib.cm import get_cmap from matplotlib.cm import get_cmap
from mpl_toolkits.basemap import Basemap from mpl_toolkits.basemap import Basemap
from wrf import to_np, getvar, smooth2d, latlon_coords, get_basemap from wrf import to_np, getvar, smooth2d, get_basemap, latlon_coords
# Open the NetCDF file
ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00") ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00")
# Get the sea level pressure # Get the sea level pressure
slp = getvar(ncfile, "slp") slp = getvar(ncfile, "slp")
# Smooth the sea level pressure since it tends to be noise y near the mountains # Smooth the sea level pressure since it tends to be noisy near the mountains
smooth_slp = smooth2d(slp, 3) smooth_slp = smooth2d(slp, 3)
# Get the latitude and longitude coordinate s # Get the latitude and longitude point s
lats, lons = latlon_coords(slp) lats, lons = latlon_coords(slp)
# Get the basemap projection object # Get the basemap object
bm = get_basemap(slp) bm = get_basemap(slp)
# Create a figure that's 10x10 # Create a figure
fig = plt.figure(figsize=(10,10 )) fig = plt.figure(figsize=(12,9 ))
# Add geographic outlines
bm.drawcoastlines(linewidth=0.25) bm.drawcoastlines(linewidth=0.25)
bm.drawstates(linewidth=0.25) bm.drawstates(linewidth=0.25)
bm.drawcountries(linewidth=0.25) bm.drawcountries(linewidth=0.25)
@ -362,12 +355,13 @@ Plotting a Two-Dimensional Field
# Draw the contours and filled contours # Draw the contours and filled contours
bm.contour(x, y, to_np(smooth_slp), 10, colors="black") bm.contour(x, y, to_np(smooth_slp), 10, colors="black")
bm.contourf(x, y, to_np(smooth_slp), 10) bm.contourf(x, y, to_np(smooth_slp), 10, cmap=get_cmap("jet") )
# Add a color bar # Add a color bar
plt.colorbar(shrink=.47 ) plt.colorbar(shrink=.62 )
plt.title("Sea Level Pressure (hPa)") plt.title("Sea Level Pressure (hPa)")
plt.show() plt.show()
@ -387,6 +381,7 @@ Horizontal Interpolation to a Pressure Level
from wrf import getvar, interplevel, to_np, get_basemap, latlon_coords from wrf import getvar, interplevel, to_np, get_basemap, latlon_coords
# Open the NetCDF file
ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00") ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00")
# Extract the pressure, geopotential height, and wind variables # Extract the pressure, geopotential height, and wind variables
@ -405,27 +400,28 @@ Horizontal Interpolation to a Pressure Level
# Get the lat/lon coordinates # Get the lat/lon coordinates
lats, lons = latlon_coords(ht_500) lats, lons = latlon_coords(ht_500)
# Get the basemap projection object # Get the basemap object
bm = get_basemap(ht_500) bm = get_basemap(ht_500)
# Create the figure # Create the figure
fig = plt.figure(figsize=(8,8 )) fig = plt.figure(figsize=(12,9 ))
ax = plt.axes([0.1,0.1,0.8,0.8] ) ax = plt.axes()
# Get the x and y coordinates # Convert the lat/lon coordinates to x/y coordinates in the projection space
x, y = bm(to_np(lons), to_np(lats)) x, y = bm(to_np(lons), to_np(lats))
# Create the 500 hPa geopotential height contours # Add the 500 hPa geopotential height contours
levels = np.arange(520., 580., 6.) levels = np.arange(520., 580., 6.)
contours = bm.contour(x, y, to_np(ht_500), levels=levels, colors="black") contours = bm.contour(x, y, to_np(ht_500), levels=levels, colors="black")
plt.clabel(contours, inline=1, fontsize=10, fmt="%i") plt.clabel(contours, inline=1, fontsize=10, fmt="%i")
# Wind S peed contours # Add the wind s peed contours
levels = [25, 30, 35, 40, 50, 60, 70, 80, 90, 100, 110, 120] 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, wspd_contours = bm.contourf(x, y, to_np(wspd_500), levels=levels,
cmap=get_cmap("rainbow")) cmap=get_cmap("rainbow"))
plt.colorbar(wspd_contours, ax=ax, orientation="horizontal", pad=.05) plt.colorbar(wspd_contours, ax=ax, orientation="horizontal", pad=.05)
# Add the geographic boundaries
bm.drawcoastlines(linewidth=0.25) bm.drawcoastlines(linewidth=0.25)
bm.drawstates(linewidth=0.25) bm.drawstates(linewidth=0.25)
bm.drawcountries(linewidth=0.25) bm.drawcountries(linewidth=0.25)
@ -438,6 +434,7 @@ Horizontal Interpolation to a Pressure Level
plt.show() plt.show()
Panel Plots from the Front Page Panel Plots from the Front Page
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@ -458,9 +455,8 @@ plot, see :ref:`cross_example`.
from wrf import getvar, to_np, vertcross, smooth2d, CoordPair, get_basemap, latlon_coords from wrf import getvar, to_np, vertcross, smooth2d, CoordPair, get_basemap, latlon_coords
# Open the output NetCDF file # Open the NetCDF file
filename = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-10-07_00_00_00" ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00")
ncfile = Dataset(filename)
# Get the WRF variables # Get the WRF variables
slp = getvar(ncfile, "slp") slp = getvar(ncfile, "slp")
@ -475,29 +471,28 @@ plot, see :ref:`cross_example`.
start_point = CoordPair(lat=26.76, lon=-80.0) start_point = CoordPair(lat=26.76, lon=-80.0)
end_point = CoordPair(lat=26.76, lon=-77.8) 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 # Compute the vertical cross-section interpolation. Also, include the lat/lon points along the cross-section i n
# in the metadata by setting latlon to True. # 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) 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) 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) dbz_cross = 10.0 * np.log10(z_cross)
# Extract the latitude and longitude coordinate array s # Get the latitude and longitude point s
lats, lons = latlon_coords(slp) lats, lons = latlon_coords(slp)
# Create a figure that will have 3 subplots. One on the left, # Create the figure that will have 3 subplots
# two vertically stacked on the right. fig = plt.figure(figsize=(10,7))
fig = plt.figure(figsize=(8,5))
ax_ctt = fig.add_subplot(1,2,1) ax_ctt = fig.add_subplot(1,2,1)
ax_wspd = fig.add_subplot(2,2,2) ax_wspd = fig.add_subplot(2,2,2)
ax_dbz = fig.add_subplot(2,2,4) ax_dbz = fig.add_subplot(2,2,4)
# Get the basemap pr ojection class # Get the basemap ob ject
bm = get_basemap(slp) bm = get_basemap(slp)
# Get the x, y values # Convert the lat/lon points in to x/y points in the projection space
x, y = bm(to_np(lons), to_np(lats)) x, y = bm(to_np(lons), to_np(lats))
# Make the pressure contours. # Make the pressure contours
contour_levels = [960, 965, 970, 975, 980, 990] contour_levels = [960, 965, 970, 975, 980, 990]
c1 = bm.contour(x, y, to_np(smooth_slp), levels=contour_levels, colors="white", c1 = bm.contour(x, y, to_np(smooth_slp), levels=contour_levels, colors="white",
zorder=3, linewidths=1.0, ax=ax_ctt) zorder=3, linewidths=1.0, ax=ax_ctt)
@ -511,12 +506,11 @@ plot, see :ref:`cross_example`.
bm.plot([point_x[0], point_x[1]], [point_y[0], point_y[1]], color="yellow", bm.plot([point_x[0], point_x[1]], [point_y[0], point_y[1]], color="yellow",
marker="o", zorder=3, ax=ax_ctt) marker="o", zorder=3, ax=ax_ctt)
# Create the label bar for cloud top temperature # Create the color bar for cloud top temperature
cb_ctt = fig.colorbar(ctt_contours, ax=ax_ctt, shrink=.68 ) cb_ctt = fig.colorbar(ctt_contours, ax=ax_ctt, shrink=.60 )
cb_ctt.ax.tick_params(labelsize=5) cb_ctt.ax.tick_params(labelsize=5)
# Draw the oceans, land, and states. Use the same colors as the cartopy # Draw the oceans, land, and states
# example
bm.drawcoastlines(linewidth=0.25, ax=ax_ctt) bm.drawcoastlines(linewidth=0.25, ax=ax_ctt)
bm.drawstates(linewidth=0.25, ax=ax_ctt) bm.drawstates(linewidth=0.25, ax=ax_ctt)
bm.drawcountries(linewidth=0.25, ax=ax_ctt) bm.drawcountries(linewidth=0.25, ax=ax_ctt)
@ -526,13 +520,12 @@ plot, see :ref:`cross_example`.
# Draw Parallels # Draw Parallels
parallels = np.arange(np.amin(lats), 30., 2.5) parallels = np.arange(np.amin(lats), 30., 2.5)
bm.drawparallels(parallels, ax=ax_ctt) bm.drawparallels(parallels, ax=ax_ctt, color="white" )
merids = np.arange(-85.0, -72.0, 2.5) merids = np.arange(-85.0, -72.0, 2.5)
bm.drawmeridians(merids, ax=ax_ctt) bm.drawmeridians(merids, ax=ax_ctt, color="white" )
# Get the x and y coordinate ranges for the cropped region around the # Crop the image to the hurricane region
# hurricane
x_start, y_start = bm(-85.0, np.amin(lats)) x_start, y_start = bm(-85.0, np.amin(lats))
x_end, y_end = bm(-72.0, 30.0) x_end, y_end = bm(-72.0, 30.0)
@ -540,14 +533,14 @@ plot, see :ref:`cross_example`.
ax_ctt.set_ylim([y_start, y_end]) ax_ctt.set_ylim([y_start, y_end])
# Make the contour plot for wspd # Make the contour plot for wspd
wspd_contours = ax_wspd.contourf(to_np(wspd_cross)) wspd_contours = ax_wspd.contourf(to_np(wspd_cross), cmap=get_cmap("jet") )
# Add the color bar # Add the color bar
cb_wspd = fig.colorbar(wspd_contours, ax=ax_wspd) cb_wspd = fig.colorbar(wspd_contours, ax=ax_wspd)
cb_wspd.ax.tick_params(labelsize=5) cb_wspd.ax.tick_params(labelsize=5)
# Make the contour plot for dbz # Make the contour plot for dbz
levels = [5 + 5*n for n in range(15)] levels = [5 + 5*n for n in range(15)]
dbz_contours = ax_dbz.contourf(to_np(dbz_cross), levels=levels) dbz_contours = ax_dbz.contourf(to_np(dbz_cross), levels=levels, cmap=get_cmap("jet") )
cb_dbz = fig.colorbar(dbz_contours, ax=ax_dbz) cb_dbz = fig.colorbar(dbz_contours, ax=ax_dbz)
cb_dbz.ax.tick_params(labelsize=5) cb_dbz.ax.tick_params(labelsize=5)
@ -573,7 +566,7 @@ plot, see :ref:`cross_example`.
ax_wspd.set_ylabel("Height (m)", fontsize=5) ax_wspd.set_ylabel("Height (m)", fontsize=5)
ax_dbz.set_ylabel("Height (m)", fontsize=5) ax_dbz.set_ylabel("Height (m)", fontsize=5)
# Add a title # Add titles
ax_ctt.set_title("Cloud Top Temperature (degC)", {"fontsize" : 7}) ax_ctt.set_title("Cloud Top Temperature (degC)", {"fontsize" : 7})
ax_wspd.set_title("Cross-Section of Wind Speed (kt)", {"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}) ax_dbz.set_title("Cross-Section of Reflectivity (dBZ)", {"fontsize" : 7})
@ -598,34 +591,34 @@ plotted using the standard matplotlib API.
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap from matplotlib.cm import get_cmap
import cartopy.crs as crs
from cartopy.feature import NaturalEarthFeature
from netCDF4 import Dataset from netCDF4 import Dataset
from wrf import to_np, getvar, CoordPair, vertcross from wrf import to_np, getvar, CoordPair, vertcross
# Open the output NetCDF file with PyNIO # Open the NetCDF file
filename = "wrfout_d01_2016-10-07_00_00_00" filename = "wrfout_d01_2016-10-07_00_00_00"
ncfile = Dataset(filename) ncfile = Dataset(filename)
# Extract pressure and model height # Extract the model height and wind speed
z = getvar(ncfile, "z") z = getvar(ncfile, "z")
wspd = getvar(ncfile, "uvmet_wspd_wdir", units="kt")[0,:] wspd = getvar(ncfile, "uvmet_wspd_wdir", units="kt")[0,:]
# Define the start point and end point for the cross section using # Create the start point and end point for the cross section
# latitude and longitude coordinates.
start_point = CoordPair(lat=26.76, lon=-80.0) start_point = CoordPair(lat=26.76, lon=-80.0)
end_point = CoordPair(lat=26.76, lon=-77.8) end_point = CoordPair(lat=26.76, lon=-77.8)
# Compute the vertical cross-section interpolation. Also, include the lat/lon # Compute the vertical cross-section interpolation. Also, include the lat/lon points along the cross-section.
# 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)
wspd_cross = vertcross(wspd, z, wrfin=ncfile, start_point=start_point,
end_point=end_point, latlon=True, meta=True)
# Create the figure # Create the figure
fig = plt.figure(figsize=(10,5 )) fig = plt.figure(figsize=(12,6 ))
ax = plt.axes([0.1,0.1,0.8,0.8] ) ax = plt.axes()
# Make the contour plot # Make the contour plot
wspd_contours = ax.contourf(to_np(wspd_cross)) wspd_contours = ax.contourf(to_np(wspd_cross), cmap=get_cmap("jet"))
# Add the color bar # Add the color bar
plt.colorbar(wspd_contours, ax=ax) plt.colorbar(wspd_contours, ax=ax)
Bill Ladwig
8 years ago
