Caribbean precipitation extremes using HighResMIP¶

Overview¶

In this notebook, we’ll do the following:

Learn how to make subregional-scale plots with the CMIP6 model output
Plot spatial trends in summertime environmental variables
Calculate the genesis potential index (GPI), an indicator of environmental favorability for tropical cyclogenesis

Prerequisites¶

Concepts	Importance	Notes
Intro to Cartopy	Necessary
Load CMIP6 data with Intake-ESM	Necessary
Intro to Xarray	Necessary
Dask arrays with Xarray	Helpful

Time to learn: 30 minutes.

Installs¶

Before you start, install the xMIP, regionmask, and tcpyPI packages. These packages allow us to mask and handle the bulky CMIP6 model output with ease. tcpyPI helps us calculate the potential intensity of storms given certain environmental conditions. Remove -q below for pip install output.

# !pip install -q git+https://github.com/jbusecke/xmip.git

# !pip install -q regionmask==0.12.1

# !pip install -q tcpypi

Imports¶

Begin your body of content with another --- divider before continuing into this section, then remove this body text and populate the following code cell with all necessary Python imports up-front:

%matplotlib inline
from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
import xarray as xr
import regionmask
import intake
import xesmf as xe
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from xmip.preprocessing import combined_preprocessing
from xmip.regionmask import merged_mask
from xgcm import Grid
import xgcm 
from xgcm.autogenerate import generate_grid_ds

from cartopy.mpl.geoaxes import GeoAxes
import matplotlib.ticker as mticker

plt.rcParams['figure.figsize'] = 12, 6
plt.rcParams['figure.titlesize'] = 15
plt.rcParams["axes.titlesize"] = 15

Initiating compute cluster¶

import os
import sys
from dask_gateway import Gateway
from dask.distributed import Client


platform = sys.platform

if (platform == 'win32'):
    import multiprocessing.popen_spawn_win32
else:
    import multiprocessing.popen_spawn_posix

client = Client()
client

Once you initiate your Client, navigate to the Dask <Dask> tab in the far-left sidebar or click the Launch dashboard in JupyterLab button above. In the search bar at the top, input your cluster’s url. Then, replace https://127.0.0.1: with /proxy/.

Load the CMIP6 HighResMIP data¶

Since we’re examining precipitation extremes at regional scales, it’s good to use the highest resolution data available. In this notebook, we’ll use the HighResMIP experiments, that provide a relatively higher resolution than the CMIP6 data. In this section, we’ll choose one global climate model (as an example, NOAA GFDL’s GFDL-CM4C192 model) for our analyses. We’ll grab the following environmental variables:

Omega (wap)
Convective precipitation (prc)
Mean sea-level pressure (psl)
Eastward, Meridional (horizontal) winds (ua, va)
Relative humidity (hur)
Specific humidity (hus)

Below, we’re using Intake-ESM to grab the model, experiments, and variables we want. Intake-ESM packages all of our data into two datasets (by experiment).

cat_url = "https://storage.googleapis.com/cmip6/pangeo-cmip6.json"
col = intake.open_esm_datastore(cat_url)
col

Use the column headers above to refine your catalog search. Here’s a handy list of CMIP6 experiments and variables to help your search.

cat = col.search(activity_id='HighResMIP', 
                 institution_id='NOAA-GFDL', 
                 experiment_id=['highresSST-present', 'highresSST-future'],
                 table_id='Amon',
                 variable_id=['psl', 'tas', 'wap', 'prc', 'hur', 'ua', 'va', 'hus'])

cat.df

dset_dict = cat.to_dataset_dict(zarr_kwargs={'consolidated': True})
list(dset_dict.keys())

We’ll consolidate our catalog into two datasets, then merge with xarray.concat to create a larger dataset.

dsPRES = dset_dict['HighResMIP.NOAA-GFDL.GFDL-CM4C192.highresSST-present.Amon.gr3'].squeeze()
dsFUT = dset_dict['HighResMIP.NOAA-GFDL.GFDL-CM4C192.highresSST-future.Amon.gr3'].squeeze()
dsFULL = xr.concat([dsPRES, dsFUT], dim="time")
dsFULL

The global warming trend and spatial trend maps¶

We understand that the Earth is warming from increasing linear trends in global mean surface temperatures. But, not everywhere warms at the same rate. This is known as the pattern effect, and it results in more complex global and regional climate change impacts. Using xarray, we can find the linear trend at each gridpoint in our data.

Let’s calculate the trend in precipitation across the Caribbean for our GCM.

import matplotlib.gridspec as gridspec
from cartopy.mpl.ticker import LatitudeFormatter, LongitudeFormatter

# Calculate the global mean time series and linear trend from surface air temperatures
regional_mean = dsFULL.prc.sel(time=dsFULL.time.dt.month.isin([8,9,10]), lat=slice(10, 30), lon=slice(260, 310)).mean(['lat', 'lon'])
regional_anom = regional_mean.groupby('time.month') - regional_mean.groupby('time.month').mean('time')
regional_anom = regional_anom*1e4

# Use numpy lstsq to find linear trend
x = np.arange(0, regional_anom.time.shape[0])
A = np.vstack([x, np.ones(len(x))]).T

# Calculate the least squares fit
m, c = np.linalg.lstsq(A, regional_anom.values, rcond=None)[0]
prc_linear = m*x + c

# Calculate the linear trend at each grid point
da = dsFULL.prc.sel(time=dsFULL.time.dt.month.isin([8,9,10]), lat=slice(10, 30), lon=slice(260, 310)).groupby('time.year').mean('time')*1e4
res = da.polyfit(dim="year", deg=1)

# Plot results
fig = plt.figure(figsize=(12,10), dpi=200)

gs = gridspec.GridSpec(2, 4, height_ratios=[0.35, 0.65], hspace=-0.01) 

ax1= fig.add_subplot(gs[0, 0:3])
ax2 = fig.add_subplot(gs[1:2, 0:3], projection=ccrs.PlateCarree())

# First subplot
regional_anom.plot(ax=ax1, linewidth=1, color='black')
ax1.plot(regional_anom.time, prc_linear, color='red')
ax1.set_title('Regional mean Precipitation over the Caribbean')

# Second subplot
ax2.set_extent([-100, -50, 10, 30], crs=ccrs.PlateCarree())
ax2.add_feature(cfeature.LAND, zorder=2, facecolor='None', edgecolor='black')

ax2.set_xticks(np.arange(-100, -50, 20), crs=ccrs.PlateCarree())
ax2.set_yticks(np.arange(10, 35, 5), crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax2.xaxis.set_major_formatter(lon_formatter)
ax2.yaxis.set_major_formatter(lat_formatter)

fc = res.polyfit_coefficients.sel(degree=1).plot(ax=ax2, cmap=plt.cm.BrBG, extend='both', add_colorbar=False)
ax2.set_title('Linear trend (i.e. slope) from 1950-2050 at each gridpoint')

cax1 = fig.add_axes([0.73, 0.235, 0.02, 0.25])
cbar = fig.colorbar(fc, cax=cax1, extend='both')

/home/runner/micromamba/envs/cookbook-dev/lib/python3.13/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/50m_physical/ne_50m_land.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

In our model, the Gulf of Mexico and western Caribbean regions experiences more drying relative to the eastern Caribbean. Another useful field to examine is the 500-hPa omega ( $\omega500$ field, that tells us where deep convection occurs. Let’s rerun the cells above to examine $\omega500$ trends. Note that, in contrast to precipitation, $\omega500$ is negative for upward motion (convective activity) and positive to downward motion (no convective activity).

# Calculate the global mean time series and linear trend from surface air temperatures
regional_mean = dsFULL.wap.sel(time=dsFULL.time.dt.month.isin([8,9,10]), lat=slice(10, 30), lon=slice(260, 310), plev=500*100).mean(['lat', 'lon'])
regional_anom = regional_mean.groupby('time.month') - regional_mean.groupby('time.month').mean('time')
regional_anom = regional_anom*1e2

# Use numpy lstsq to find linear trend
x = np.arange(0, regional_anom.time.shape[0])
A = np.vstack([x, np.ones(len(x))]).T

# Calculate the least squares fit
m, c = np.linalg.lstsq(A, regional_anom.values, rcond=None)[0]
prc_linear = m*x + c

# Calculate the linear trend at each grid point
da = dsFULL.wap.sel(time=dsFULL.time.dt.month.isin([8,9,10]), lat=slice(10, 30), lon=slice(260, 310), plev=500*100).groupby('time.year').mean('time')*1e2
res = da.polyfit(dim="year", deg=1)

# Plot results
fig = plt.figure(figsize=(12,10), dpi=200)

gs = gridspec.GridSpec(2, 4, height_ratios=[0.35, 0.65], hspace=-0.02) 

ax1= fig.add_subplot(gs[0, 0:3])
ax2 = fig.add_subplot(gs[1:2, 0:3], projection=ccrs.PlateCarree())

# First subplot
regional_anom.plot(ax=ax1, linewidth=1, color='black')
ax1.plot(regional_anom.time, prc_linear, color='red')
ax1.set_title('Regional mean Omega over the Caribbean')

# Second subplot
ax2.set_extent([-100, -50, 10, 30], crs=ccrs.PlateCarree())
ax2.add_feature(cfeature.LAND, zorder=2, facecolor='None', edgecolor='black')

ax2.set_xticks(np.arange(-100, -50, 20), crs=ccrs.PlateCarree())
ax2.set_yticks(np.arange(10, 35, 5), crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax2.xaxis.set_major_formatter(lon_formatter)
ax2.yaxis.set_major_formatter(lat_formatter)

fc = res.polyfit_coefficients.sel(degree=1).plot(ax=ax2, cmap=plt.cm.RdBu_r, extend='both', add_colorbar=False)
ax2.set_title('Linear trend (i.e. slope) from 1950-2050 at each gridpoint')

cax1 = fig.add_axes([0.73, 0.235, 0.02, 0.25])
cbar = fig.colorbar(fc, cax=cax1, extend='both')

Environmental changes by experiment¶

Here, we look at the simulated 500-hPa omega ( $\omega500$ ) fields for the highresSST-Present and highresSST-Future. The $\omega500$ variable is a useful indicator for examining where extreme convection tends to occur.

# Subset our data to only look at months August-October
da1_env = dsPRES.wap.sel(time=dsPRES.time.dt.month.isin([8,9,10]), plev=500*100).mean('time')*1e2
da2_env = dsFUT.wap.sel(time=dsFUT.time.dt.month.isin([8,9,10]), plev=500*100).mean('time')*1e2

difference = da2_env - da1_env

fig = plt.figure(figsize=(12,10), dpi=200)

ax = fig.add_subplot(3, 1, 1, projection=ccrs.PlateCarree())
ax.set_extent([-100, -50, 10, 30], crs=ccrs.PlateCarree())
ax.add_feature(cfeature.LAND, zorder=2, facecolor='None', edgecolor='black')

ax.set_xticks(np.arange(-100, -50, 20), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(10, 40, 10), crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)

fc = da1_env.plot(levels=np.arange(-5, 5.5, 0.5), cmap=plt.cm.RdGy_r, extend='both', add_colorbar=False, ax=ax)
ax.set_title(r"$\omega500$; highresSST-Present: 1950-2014")

cbar = plt.colorbar(fc, orientation='vertical', ax=ax)
cbar.ax.tick_params(labelsize=16)
cbar.ax.set_title(r'Pa s$^{-1}$')

ax = fig.add_subplot(3, 1, 2, projection=ccrs.PlateCarree())
ax.set_extent([-100, -50, 10, 30], crs=ccrs.PlateCarree())
ax.add_feature(cfeature.LAND, zorder=2, facecolor='None', edgecolor='black')

ax.set_xticks(np.arange(-100, -50, 20), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(10, 40, 10), crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)

fc = da2_env.plot(levels=np.arange(-5, 5.5, 0.5), cmap=plt.cm.RdGy_r, extend='both', add_colorbar=False, ax=ax)
ax.set_title(r"$\omega500$:highresSST-Future: 2015-2050")

cbar = plt.colorbar(fc, orientation='vertical', ax=ax)
cbar.ax.tick_params(labelsize=16)
cbar.ax.set_title(r'Pa s$^{-1}$')
        
ax = fig.add_subplot(3, 1, 3, projection=ccrs.PlateCarree())
ax.set_extent([-100, -50, 10, 30], crs=ccrs.PlateCarree())
ax.add_feature(cfeature.LAND, zorder=2, facecolor='None', edgecolor='black')

ax.set_xticks(np.arange(-100, -50, 20), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(10, 40, 10), crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)

fc = difference.plot(levels=np.arange(-3, 3.2, 0.2), cmap=plt.cm.RdBu_r, extend='both', add_colorbar=False, ax=ax)
ax.set_title(r"$\omega500$:highresSST-Future - highresSST-Present")

cbar = plt.colorbar(fc, orientation='vertical', ax=ax)
cbar.ax.tick_params(labelsize=16)
cbar.ax.set_title(r'Pa s$^{-1}$')

plt.subplots_adjust(hspace=0.4)

Calculating genesis potential¶

Here, we’ll do something a little different. We’ll retrieve some more data above to calculate what is known as genesis potential for the Caribbean region. Genesis potential is an indicator of how favorable

client.close()

Summary¶

In the above workflow, we examine environmental changes in the CMIP6 HighResMIP simulations and calculate changes in GPI, an indicator of tropical cyclone activity for the Caribbean region.

What’s next?¶

In the next book, we examine how El Niño Southern Oscillation (ENSO) variability is a key factor in this regional-scale variability and trends.

Resources and references¶

Finally, be rigorous in your citations and references as necessary. Give credit where credit is due. Also, feel free to link to relevant external material, further reading, documentation, etc. Then you’re done! Give yourself a quick review, a high five, and send us a pull request. A few final notes:

Kernel > Restart Kernel and Run All Cells... to confirm that your notebook will cleanly run from start to finish
Kernel > Restart Kernel and Clear All Outputs... before committing your notebook, our machines will do the heavy lifting
Take credit! Provide author contact information if you’d like; if so, consider adding information here at the bottom of your notebook
Give credit! Attribute appropriate authorship for referenced code, information, images, etc.
Only include what you’re legally allowed: no copyright infringement or plagiarism

Thank you for your contribution!

Notebook Gallery

Extreme SSTs