Contour Maps of Hurricane Maria on a 15–3 km MPAS Mesh

Contour Maps of Hurricane Maria on a 15–3 km MPAS Mesh#

Author: Dylan Nelson

This example draws contour maps from a regional MPAS-A simulation of September 2017, a few hours before Hurricane Maria made landfall in Puerto Rico:

  1. Sea-level pressure contours over 10 m wind speed around the storm.

  2. 500 hPa geopotential height contours over precipitable water across the region.

The contour lines are drawn directly on the unstructured mesh, without regridding, with .plot.contour(). Contour Demos has more examples of the method.

The last three sections go beyond maps and use other UXarray functions on the same data: an azimuthal mean, a cross-section and a weighted mean.

Imports#

import warnings

import cartopy.crs as ccrs
import geoviews.feature as gf
import holoviews as hv
import numpy as np

import uxarray as ux

warnings.filterwarnings("ignore", message=r".*(No spatial coordinate|np\.dot\(\))")

Data#

The simulation is described in Regional Multi-Week Convection-Permitting Simulations using the Model for Prediction Across Scales (MPAS-A) for use in the Caribbean, Mexico, and Central America Regions (Núñez Ocasio, Bacmeister, Xue and Morrison; NSF NCAR GDEX dataset d010077, doi:10.5065/PNB8-GD97)

grid_path = "/gdex/data/d010077/native_diag/MAAG.static.nc"
data_path = "/gdex/data/d010077/native_diag/diag.2017-09-20_06.00.00.nc"

uxds = ux.open_dataset(grid_path, data_path)

Preparing the Fields#

Computing and transforming the variables we want to use

mslp = uxds["mslp"].isel(Time=0) / 100  # sea-level pressure, Pa converted to hPa
u10 = uxds["u10"].isel(Time=0)  # 10 m wind, eastward component (m/s)
v10 = uxds["v10"].isel(Time=0)  # 10 m wind, northward component (m/s)
wind = ((u10**2 + v10**2) ** 0.5).rename("wind")  # 10 m wind speed (m/s)
z500 = uxds["height_500hPa"].isel(Time=0)  # height of the 500 hPa pressure surface (m)
pw = uxds["precipw"].isel(Time=0)  # precipitable water (kg/m², the same as mm)

Map 1: Pressure over Wind Speed#

Selecting the Storm#

bounding_box() keeps only the cells in a latitude–longitude box around the storm. Reducing the data size by subsetting will speed up the following computation steps.

storm_box = dict(lon_bounds=(-71, -59), lat_bounds=(13, 22))

mslp_storm = mslp.subset.bounding_box(**storm_box)
wind_storm = wind.subset.bounding_box(**storm_box)
print(
    f"{mslp.uxgrid.n_face:,} cells in the mesh, {mslp_storm.uxgrid.n_face:,} around the storm"
)

6,087,575 cells in the mesh, 152,236 around the storm

The map is built from three parts, each drawn alone first and then combined.

Part A - Shading with .plot()#

plot() shades the mesh cells with the wind speed. The plot is interactive: the toolbar zooms and pans, and pointing at the map shows the wind speed.

The map covers 10° of longitude and 7° of latitude, so the plot is 800 by 560 pixels. With the same proportions as the map, the view ends where the shading ends.

projection = ccrs.PlateCarree()

shading = wind_storm.plot(
    projection=projection,
    cmap="YlOrRd",
    clim=(0, 70),
    clabel="10 m wind speed (m/s)",
    xlim=(-70, -60),
    ylim=(14, 21),
    width=800,
    height=560,
    title="Part A - 10 m wind speed",
)

shading

Part B - Contour Lines with .plot.contour()#

contour() draws the pressure contours. The lines are computed on the mesh itself, so nothing has to be regridded: they follow the edges of the cells that separate pressures above a level from pressures below it. Each line is labeled with its level. method="interpolated" gives smooth lines instead.

mslp_levels = np.arange(900, 1016, 4)

contours = mslp_storm.plot.contour(
    levels=mslp_levels,
    projection=projection,
    color="black",
    line_width=1,
    xlim=(-70, -60),
    ylim=(14, 21),
    width=800,
    height=560,
    title="Part B - Sea-level pressure (hPa)",
)

contours

Part C - Coastlines with GeoViews#

The coastlines come from GeoViews. scale="50m" selects data detailed enough for the islands, and the projection and extent are the same as in the other two parts.

coastline = gf.coastline(
    scale="50m",
    projection=projection,
    xlim=(-70, -60),
    ylim=(14, 21),
    width=800,
    height=560,
    title="Part C - Coastlines",
)

coastline

Combining A, B, and C#

* overlays the three parts. They share the map projection, so they line up. Pointing at a contour line shows its pressure. It’s as easy as multiplying the variables from the past 3 plots

combined = shading * contours * coastline

combined.opts(
    title="Combining plots A, B, and C",
)

Adding Styling#

The plots so far use the default look. .opts changes how a plot is drawn after it has been made, without computing anything again. Options are grouped by the kind of layer they apply to:

  • hv.opts.Contours for the contour lines, here a slightly thicker line.

  • hv.opts.Labels for their labels, here white, bold and larger, with a black outline that keeps them readable on light and dark shading.

  • hv.opts.Feature for the coastlines, here a thicker, green line.

  • hv.opts.Overlay for the plot as a whole, here the title and its size. backend_opts reaches properties of the Bokeh plot that have no option of their own, such as the font of the title.

The number of lines is not a style. levels sets it when the contours are computed, so .plot.contour() is called again, with seven chosen levels between 950 and 1013 hPa instead of one every 4 hPa. They are wider apart near the center of the storm, where the pressure falls quickly, and closer together further out. That leaves room for the labels.

Applying Customizations in the HoloViews documentation describes .opts and how to list the options an element has.

fewer_contours = mslp_storm.plot.contour(
    levels=np.array([950, 970, 990, 1000, 1005, 1010, 1013]),
    projection=projection,
    color="black",
)

styled = shading * fewer_contours * coastline

styled.opts(
    hv.opts.Contours(
        line_width=1.2,
        line_color="black",
    ),
    hv.opts.Labels(
        text_color="white",
        text_font_size="16pt",
        text_font_style="bold",
        text_outline_color="black",
    ),
    hv.opts.Feature(
        line_width=1,
        line_color="green",
    ),
    hv.opts.Overlay(
        title="Hurricane Maria, 06 UTC 20 September 2017\nSea-level pressure (hPa) over 10 m wind speed",
        fontsize={"title": "16pt"},
        backend_opts={"plot.title.text_font": "Georgia"},
    ),
)

Map 2: 500 hPa Height over Precipitable Water#

The same three parts at a larger scale, covering most of the simulation domain, with a different colormap and contour levels.

region_box = dict(lon_bounds=(-105, -35), lat_bounds=(0, 48))

z500_region = z500.subset.bounding_box(**region_box)
pw_region = pw.subset.bounding_box(**region_box)

shading = pw_region.plot(
    projection=projection,
    cmap="GnBu",
    clim=(10, 70),
    clabel="Precipitable water (mm)",
    xlim=(-100, -40),
    ylim=(5, 40),
    width=900,
    height=525,
    title="06 UTC 20 September 2017: 500 hPa geopotential height (m) and precipitable water",
)
contours = z500_region.plot.contour(
    levels=np.arange(5700, 5950, 30),
    projection=projection,
    color="black",
    line_width=1,
)

shading * contours * gf.coastline(scale="50m", projection=projection)

Plot 3: Azimuthal Mean around the Eye#

The eye of the storm is taken to be the cell with the lowest sea-level pressure.

eye = int(mslp_storm.argmin(dim="n_face"))
eye_lon = float(mslp_storm.uxgrid.face_lon[eye])
eye_lat = float(mslp_storm.uxgrid.face_lat[eye])

print(
    f"Lowest pressure: {float(mslp_storm.min()):.0f} hPa at {eye_lat:.2f}°N, {abs(eye_lon):.2f}°W"
)

Lowest pressure: 948 hPa at 16.21°N, 65.23°W

azimuthal_mean() averages a field along circles around a point, here the eye, every 0.05° out to 3° (1° of arc is about 111 km). The result is an ordinary Xarray array with a radius dimension. It is drawn as two HoloViews curves: wind speed and pressure against the distance from the eye.

wind_profile = wind_storm.azimuthal_mean(
    center_coord=(eye_lon, eye_lat),
    outer_radius=3,
    radius_step=0.05,
)
mslp_profile = mslp_storm.azimuthal_mean(
    center_coord=(eye_lon, eye_lat),
    outer_radius=3,
    radius_step=0.05,
)

distance = hv.Dimension("radius", label="Distance from the eye (degrees)")

wind_curve = hv.Curve(
    (wind_profile.radius, wind_profile),
    distance,
    hv.Dimension("wind", label="10 m wind speed (m/s)"),
).opts(
    title="Plot 3 - Wind speed around the eye",
    width=400,
    height=320,
)
mslp_curve = hv.Curve(
    (mslp_profile.radius, mslp_profile),
    distance,
    hv.Dimension("mslp", label="Sea-level pressure (hPa)"),
).opts(
    title="Plot 3 - Pressure around the eye",
    width=400,
    height=320,
)

wind_curve + mslp_curve

Plot 4: Zonal Cross-Section through the Eye#

cross_section() samples a field along a line, here the line of constant latitude through the eye. With lat= the samples cover the whole latitude circle, so steps=14401 places one every 0.025° of longitude, about the width of a mesh cell, and dropna keeps the samples that fall on the cells around the storm.

wind_section = wind_storm.cross_section(
    lat=eye_lat,
    steps=14401,
).dropna("steps")
mslp_section = mslp_storm.cross_section(
    lat=eye_lat,
    steps=14401,
).dropna("steps")

longitude = hv.Dimension("lon", label="Longitude")

wind_curve = hv.Curve(
    (wind_section.lon, wind_section),
    longitude,
    hv.Dimension("wind", label="10 m wind speed (m/s)"),
).opts(
    title=f"Plot 4 - Wind speed along {eye_lat:.2f}°N",
    xlim=(-70, -60),
    width=400,
    height=320,
)
mslp_curve = hv.Curve(
    (mslp_section.lon, mslp_section),
    longitude,
    hv.Dimension("mslp", label="Sea-level pressure (hPa)"),
).opts(
    title=f"Plot 4 - Pressure along {eye_lat:.2f}°N",
    xlim=(-70, -60),
    width=400,
    height=320,
)

wind_curve + mslp_curve