forked from 3rdparty/wrf-python
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.
1.2 MiB
1.2 MiB
Cartopy Examples¶
In [1]:
%matplotlib inline
In [2]:
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 to_np, getvar, smooth2d, get_cartopy, latlon_coords
# 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
lats, lons = latlon_coords(slp)
# Get the cartopy projection class
cart_proj = get_cartopy(slp)
# Create the figure
fig = plt.figure(figsize=(8,8))
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, to_np(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]
contour_levels = np.arange(-80.0, 10.0, 10.0)
plt.contourf(to_np(lons), to_np(lats), to_np(ctt), contour_levels, cmap=get_cmap("Greys"),
transform=crs.PlateCarree(), zorder=2)
plt.plot([-80,-77.8], [26.76,26.76], 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, shrink=.8)
# 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 [3]:
%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 to_np, getvar, smooth2d, ll_to_xy, CoordPair, vertcross, getproj, get_proj_params, to_xy_coords, latlon_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(to_np(wspd_cross))
# Add the color bar
fig.colorbar(a, ax=axes[0])
b = axes[1].contourf(to_np(dbz_cross), levels=levels)
fig.colorbar(b, ax=axes[1])
# 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)]
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 = to_np(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 [4]:
%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, to_np, vertcross, smooth2d, CoordPair, get_cartopy, latlon_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)
# 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 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
lats, lons = latlon_coords(slp)
# Get the cartopy projection class
cart_proj = get_cartopy(slp)
# 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, 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(to_np(lons), to_np(lats), to_np(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(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.savefig("/Users/ladwig/Documents/workspace/wrf_python/doc/source/_static/images/matthew.png",
transparent=True, bbox_inches="tight")
plt.show()
In [5]:
%matplotlib inline
# 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 to_np, getvar, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim
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, lons = latlon_coords(slp)
# 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(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. Not really necessary, but used for demonstration.
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)")
Out[5]:
In [6]:
%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 to_np, getvar, CoordPair, vertcross
# 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")
wspd = getvar(ncfile, "uvmet_wspd_wdir", units="kt")[0,:]
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.
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()
In [7]:
from __future__ import (absolute_import, division, print_function, unicode_literals)
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, latlon_coords, get_cartopy, cartopy_xlim, cartopy_ylim
ncfile = Dataset("/Users/ladwig/Documents/wrf_files/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 map projection information
cart_proj = get_cartopy(ht_500)
# Create the figure
fig = plt.figure(figsize=(8,8))
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()
Basemap Examples¶
In [8]:
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
from mpl_toolkits.basemap import Basemap
from wrf import to_np, getvar, smooth2d, get_basemap, latlon_coords
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, lons = latlon_coords(slp)
# The basemap method returns a Basemap object
bm = get_basemap(slp)
# Create a figure that's 8x8
fig = plt.figure(figsize=(8,8))
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()
In [9]:
%matplotlib inline
from __future__ import (absolute_import, division, print_function, unicode_literals)
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("/Users/ladwig/Documents/wrf_files/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 map projection information
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()
In [10]:
%matplotlib inline
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 with PyNIO
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 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\\
lats, lons = latlon_coords(slp)
# Create the figure which will have 3 subplots
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 cartopy 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=.75)
cb_ctt.ax.tick_params(labelsize=5)
# Draw the oceans, land, and states
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)-1.0, 30., 2.5)
bm.drawparallels(parallels, ax=ax_ctt)
merids = np.arange(-85.0, -72.0, 2.5)
bm.drawmeridians(merids, ax=ax_ctt)
# Crop the image
x_start, y_start = bm(-85.0, np.amin(lats)-1.0)
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])
## Plot the cross sections
# 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.savefig("/Users/ladwig/Documents/workspace/wrf_python/doc/source/_static/images/matthew.png",
# transparent=True, bbox_inches="tight")
plt.show()
In [11]:
# SLP
from __future__ import (absolute_import, division, print_function, unicode_literals)
import Ngl
import Nio
import numpy as np
from wrf import to_np, getvar, smooth2d, get_pyngl, latlon_coords
ncfile = Nio.open_file(b"/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-10-07_00_00_00.nc")
# Get the sea level pressure
ctt = getvar(ncfile, "ctt")
wks = Ngl.open_wks(b"png", b"test")
# Get the map projection
resources = get_pyngl(ctt)
# This needs to be False if you want to crop
resources.tfDoNDCOverlay = False
resources.mpLeftCornerLatF = 30.0
resources.mpRightCornerLatF = np.amin(to_np(ctt.coords["XLAT"]))
resources.mpLeftCornerLonF = -85.
resources.mpRightCornerLonF = -75.
lats, lons = latlon_coords(ctt, as_np=True)
resources.sfYArray = lats
resources.sfXArray = lons
resources.cnLevelSelectionMode = b"ExplicitLevels" # Define your own
resources.cnLevels = np.arange(-80.,30.,10.)
resources.cnFillOn = True
resources.cnFillMode = b"RasterFill"
#resources.cnFillPalette = Ngl.read_colormap_file(b"MPL_Greys")
Ngl.contour_map(wks, to_np(ctt), resources)
Ngl.end()
In [ ]: