Plot MPAS data using Delauney Triangulation

%matplotlib inline
import xarray as xr
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.tri as tri
import matplotlib.pyplot as plt

Get Data from Server if not locally available yet

from ngallery_utils import DATASETS
DATASETS.registry_files

%%time
filePath = DATASETS.fetch("Oklahoma.static.nc")
filePath

Read data from MPAS file

ds = xr.open_dataset(filePath)

# Remove singleton dimensions, such as Time.
ds = ds.squeeze()
ds

Convert to degrees from radians

lonData = np.degrees(ds.lonCell)
latData = np.degrees(ds.latCell)

# convert lonData to range [-180, 180]
lonData = ((lonData + 180) % 360) - 180

(np.min(lonData), np.max(lonData))

(np.min(latData), np.max(latData))

triang = tri.Triangulation(lonData, latData)

plot_data = ds["ter"]

Load colormap values and add alpha channel

from matplotlib.colors import ListedColormap
cmap_file = "./MPL_terrain.rgb"
cm_custom_values = np.loadtxt(cmap_file)
(ncolors, nchannels) = cm_custom_values.shape
(ncolors, nchannels)

cm_custom_rgba = np.hstack( (cm_custom_values, np.ones( (ncolors,1) )) )
cm_custom_rgba.shape

cmap_custom = ListedColormap(cm_custom_rgba)

Produce Figure

NOTE : First-time plotting may be slow, as in several minutes, since map features must be downloaded for the entire Earth.

fig = plt.figure(figsize=(30,30))

proj = ccrs.Orthographic(central_longitude=-97.4,central_latitude=35.1)

ax = fig.add_subplot(1,1,1, projection = proj)

# Set lat/lon bounding box and feature resolutions.
ax.set_extent([-107, -87, 25, 45])

# Choose resolution of map features.
# Note that these features are downloaded when plotting for the first time, and for the entire globe,
#  so requesting high resolution can take several minutes.
scale = '110m' # '50m' # '10m'


ax.add_feature(cfeature.LAND.with_scale(scale))
ax.add_feature(cfeature.OCEAN.with_scale(scale))
ax.add_feature(cfeature.STATES.with_scale(scale))
ax.add_feature(cfeature.LAKES.with_scale(scale), alpha=0.5)
ax.add_feature(cfeature.COASTLINE.with_scale(scale))

# Specify data range for colormap
(colormin, colormax) = (-31, 4000)

mm = ax.tripcolor(triang, plot_data, edgecolors='none', transform=ccrs.PlateCarree(), cmap=cmap_custom, vmin=colormin, vmax=colormax)

plt.colorbar(mm, orientation='horizontal', pad=0.03)
plt.title(f"MPAS terrain height ({len(ds.lonCell)} cells)", fontweight="bold", fontsize=14)

plt.show()

Save Figure to PNG File

fig.savefig('./plot_terrain.png')