NCL_WRF_zoom_1_2.py

This script illustrates the following concepts:
  • Plotting WRF data on native grid

  • Subsetting data to ‘zoom in’ on an area

  • Plotting data using wrf python functions

  • Following best practices when choosing a colormap

See following URLs to see the reproduced NCL plot & script:

Import packages

import cartopy.crs as ccrs
import geocat.datafiles as gdf
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import numpy as np
from cartopy.feature import NaturalEarthFeature
from netCDF4 import Dataset
from wrf import get_cartopy, getvar, latlon_coords, to_np

Read in the data

wrfin = Dataset(
    gdf.get("netcdf_files/wrfout_d03_2012-04-22_23_00_00_subset.nc"))
td2 = getvar(wrfin, "td2")

# Set attributes for creating plot titles later
we = getattr(wrfin, 'WEST-EAST_GRID_DIMENSION')
sn = getattr(wrfin, 'SOUTH-NORTH_GRID_DIMENSION')
lvl = getattr(wrfin, 'BOTTOM-TOP_GRID_DIMENSION')
dis = getattr(wrfin, 'DY') / 1000  # Divide by 1000 to go from m to km
phys = getattr(wrfin, 'MP_PHYSICS')
pbl = getattr(wrfin, 'BL_PBL_PHYSICS')
cu = getattr(wrfin, 'CU_PHYSICS')
s_date = getattr(wrfin, 'START_DATE')
str_format = "WE={}; SN={}; Levels={}; Dis={}km; Phys Opt={}; PBL Opt={}; Cu Opt={}"
sd_frmt = "Init: {}"

Create a subset of the data for zoomed in projection

dims = td2.shape

y_start = int(dims[0] / 2)
y_end = int(dims[0] - 1)

x_start = int(0)
x_end = int(dims[1] / 2)

td2_zoom = td2[y_start:y_end, x_start:x_end]

# Define the latitude and longitude coordinates
lats, lons = latlon_coords(td2_zoom)

Plot the data

# The `get_cartopy` wrf function will automatically find and use the
# intended map projection for this dataset
cart_proj = get_cartopy(td2_zoom)

fig = plt.figure(figsize=(12, 12))
ax = plt.axes(projection=cart_proj)

# Add features to the projection
states = NaturalEarthFeature(category="cultural",
                             scale="50m",
                             facecolor="none",
                             name="admin_1_states_provinces")

ax.add_feature(states, linewidth=0.5, edgecolor="black")
ax.coastlines('50m', linewidth=0.8)

# Add filled dew point temperature contours
plt.contourf(to_np(lons),
             to_np(lats),
             to_np(td2_zoom),
             levels=13,
             cmap="magma",
             transform=ccrs.PlateCarree(),
             vmin=-8,
             vmax=18)

# Add a colorbar
cbar = plt.colorbar(ax=ax,
                    orientation="horizontal",
                    ticks=np.arange(-6, 18, 2),
                    drawedges=True,
                    extendrect=True,
                    pad=0.08,
                    shrink=0.75,
                    aspect=30)

# Format location of colorbar text to look like NCL version
cbar.ax.text(0.5,
             1.5,
             '2m Dewpoint Temperature (C)',
             horizontalalignment='center',
             verticalalignment='center',
             transform=cbar.ax.transAxes)

# Format colorbar ticks and labels
cbar.ax.tick_params(labelsize=10)
cbar.ax.get_xaxis().labelpad = -48

# Draw gridlines
gl = ax.gridlines(crs=ccrs.PlateCarree(),
                  draw_labels=True,
                  dms=False,
                  x_inline=False,
                  y_inline=False,
                  linewidth=1,
                  color="k",
                  alpha=0.25)

# Manipulate latitude and longitude gridline numbers and spacing
gl.top_labels = False
gl.right_labels = False
gl.xlocator = mticker.FixedLocator([-120, -121, -122, -123, -124])
gl.ylocator = mticker.FixedLocator([38, 39, 40, 41, 42])
gl.xlabel_style = {"rotation": 0, "size": 15}
gl.ylabel_style = {"rotation": 0, "size": 15}
gl.xlines = False
gl.ylines = False

# Add titles to the plot
plt.title("Zoomed in plot", loc='center', x=.13, y=1.1, size=15)
plt.title("2m Dewpoint Temperature (C)", loc='left', y=1.02, size=10)
plt.title(sd_frmt.format(s_date), loc='right', y=1.1, size=10)

# Add lower text using attributes from the dataset
fig.text(0.25, 0.1, getattr(wrfin, 'TITLE'), size=12)
fig.text(0.252,
         0.08,
         str_format.format(we, sn, lvl, dis, phys, pbl, cu),
         size=12)

plt.show()
2m Dewpoint Temperature (C), Zoomed in plot, Init: 2012-04-20_00:00:00

Total running time of the script: ( 0 minutes 3.997 seconds)

Gallery generated by Sphinx-Gallery