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(). The Contour Demos example in the UXarray documentation gallery 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)

%%time 
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)

CPU times: user 2.07 s, sys: 1min 15s, total: 1min 17s
Wall time: 1min 20s

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

.subset.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.

%%time
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
CPU times: user 31.6 s, sys: 525 ms, total: 32.2 s
Wall time: 7.73 s

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.

%%time
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

CPU times: user 5.03 s, sys: 309 ms, total: 5.34 s
Wall time: 5.95 s

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

.plot.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.

%%time
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

CPU times: user 104 ms, sys: 0 ns, total: 104 ms
Wall time: 103 ms

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.

%%time
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"},
    ),
)

CPU times: user 68.3 ms, sys: 1.71 ms, total: 70 ms
Wall time: 68.9 ms

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)

%%time
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)

CPU times: user 8min 36s, sys: 3.61 s, total: 8min 39s
Wall time: 8min 40s

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.

%%time
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

CPU times: user 902 ms, sys: 516 μs, total: 903 ms
Wall time: 902 ms

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