Morphological image processing for monitoring ocean temperature extremes

---

Hillary Scannell¹, Ryan Abernathey¹, Julius Busecke¹, David John Gagne²,
LuAnne Thompson³, and Daniel Whitt⁴

¹Lamont-Doherty Earth Observatory, Columbia University, Palisades, NY, USA

²National Center for Atmospheric Research, Boulder, CO, USA

³School of Oceanography, University of Washington, Seattle, WA, USA

⁴NASA Ames Research Center, Mountain View, CA, USA

Abstract

---

Dangerous hot-water events, called marine heatwaves, cause prolonged periods of thermal stress in the marine environment that can lead to widespread coral bleaching, harmful algal blooms, unproductive fisheries, and even economic loss. Anticipating the paths of destructive marine heatwaves remains a challenge owing to the complex spatiotemporal evolution of these events. We present a novel open source package called Ocetrac that implements morphological image processing and tracking to aid in the analysis of marine heatwave detection and movement. We discuss how this approach can be applied to monitor other extreme ocean conditions, such as deoxygenation and acidification.

Many species can cope with natural temperature fluctuations

As temperatures rise, heat extremes become more deadly

We define these hot water events as Marine Heatwaves

Marine heatwaves occur throughout the global ocean

%%HTML 
<video controls autoplay width="100%" height="100%" stype="float">
  <source src="./video/sst_pos.mp4" type="video/mp4">
</video>

Motivation

---

• Marine heatwaves don't stay in one place.


• They have complex spatial connectivity and temporal behavior.


• Local analyses may not completely characterize the evolution of these events.


• Can we build a tool to identify and track marine heatwaves?

Goals of Ocetrac

---

1. Identify marine heatwaves as 2D objects from sea surface temperature anomaly images.


2. Track marine heatwave objects in both space and time.


3. Create a new labeled dataset of marine heatwave events to better understand the evolution and patterns of marine heatwaves globally.

NOAA Optimum Interpolation Sea Surface Temperature (OISST) v2.1

monthly means from September 1981–present

0.25º longitude x 0.25º latitude grid

Data Preprocessing

Compute Sea Surface Temperature Anomalies

Extract only the extreme positive values that exceed some threshold
(e.g., 90^th percentile)

Analysis-ready, Cloud-optimized Data

• NOAA OISST v2.1 dataset is available on S3 compatible object storage
• Pangeo Forge transformed netCDF source files from NOAA into a zarr store
• Read more about Pangeo Forge: https://pangeo-forge.readthedocs.io/
• Want your dataset in the cloud? Submit a staged recipe issue: https://github.com/pangeo-forge/staged-recipes/issues

import s3fs
import xarray as xr

endpoint_url = 'https://ncsa.osn.xsede.org'
fs_osn = s3fs.S3FileSystem(anon=True, client_kwargs={'endpoint_url': endpoint_url},) 

path = "Pangeo/pangeo-forge/noaa_oisst/v2.1-avhrr.zarr"
ds = xr.open_zarr(fs_osn.get_mapper(path), consolidated=True).resample(time='MS').mean()

Create and Connect to Dask Distributed Cluster

from dask_gateway import Gateway
from dask.distributed import Client

gateway = Gateway()
cluster = gateway.new_cluster() # new Dask cluster
cluster.adapt(minimum=1, maximum=30) # cluster will automatically resize itself based on the workload
cluster

client = Client(cluster) # creates a client so computations using Dask will be executed on the cluster.
client

Client Scheduler: gateway://traefik-gcp-uscentral1b-prod-dask-gateway.prod:80/prod.70ffdf5b78924f028e89ffcf15a254f1 Dashboard: /services/dask-gateway/clusters/prod.70ffdf5b78924f028e89ffcf15a254f1/status

Cluster Workers: 0 Cores: 0 Memory: 0 B

Steps in Data Preprocessing:

• Calculate the sea surface temperature (SST) anomaly as the daily SST minus the 1971-2020 climatological mean.
• Define a threshold to identify marine heatwaves. We will use the 90^th percentile computed across the entire dataset (1981-2021).
• Find where the SST anomalies exceed this threshold.
• Convert the final query into a binary image for Ocetrac.

Find where anomalies exceed the 90^th percentile threshold

import numpy as np

if ds.anom.chunks:
    ds['anom'] = ds.anom.chunk({'time': -1})
    
anom_threshold = ds.anom.quantile(.90, dim=('time'))
anom_mhw = ds.anom.where(ds.anom>=anom_threshold, other=np.nan).isel(zlev=0)

Convert the marine heatwave anomalies into a binary image

mhw_binary = anom_mhw.where(anom_mhw>0, other=0)
mhw_binary = mhw_binary.where(mhw_binary==0, other=1)

How does Ocetrac detect objects?

Multidimensional Image Processing [scipy.ndimage]

Once objects are identified, they are sorted based on area

Objects below a certain threshold are eliminated

Multiple Object Tracking

Marine heatwaves are connected when they share a surface in either x, y, or t

Marine heatwaves may split or merge

---

Ocetrac Overview

Run Ocetrac.Tracker

# Set model parameters
da = mhw_binary[:100,:,:].load() # taking a slice of the dataset
radius = 8 # radius for structuring element
min_size_quartile = 0.75 # threshold for object areas
xdim = 'lon'
ydim = 'lat'

import ocetrac.model as oce

Tracker =  oce.Tracker(da, mask, radius, min_size_quartile, xdim, ydim)
blobs = Tracker.track()

minimum area:  1901.0
inital objects identified 	 1901
final objects tracked 	 194

Resulting Marine Heatwave Objects

import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
import cartopy.crs as ccrs
import cartopy.feature as cfeature

maxl = int(np.nanmax(blobs.values))
cm = ListedColormap(np.random.random(size=(maxl, 3)).tolist())

plt.rc('font', size=20); plt.figure(figsize=(15,6)); 
ax = plt.axes(projection=ccrs.PlateCarree())
blobs.isel(time=0).plot(transform=ccrs.PlateCarree(), cmap=cm, add_colorbar=False)
ax.coastlines(resolution='110m', color='black', linewidth=1)
ax.add_feature(cfeature.LAND, facecolor='w'); 

Example marine heatwave evolution in the eastern Indian Ocean

id = 35
event = blobs.where(blobs==id, drop=True)

plt.rc('font', size=12); plt.figure(figsize=(15,6));
for i in enumerate(range(1,len(event.time))):
    ax = plt.subplot(2,3,i[1],projection=ccrs.PlateCarree())
    event.isel(time=i[0]).plot(transform=ccrs.PlateCarree(), 
                               cmap=cm, add_colorbar=False, add_labels=False)
    plt.title(event.isel(time=i[0]).time.values.astype('datetime64[D]'))
    ax.coastlines(resolution='110m', color='black', linewidth=1) 
    ax.add_feature(cfeature.LAND, facecolor='w');

# Create intensity image for the Indian Ocean marine heatwave
event_intensity = ds.anom.isel(zlev=0).where((ds.time==event.time) & 
                                             (ds.lat==event.lat) & 
                                             (ds.lon==event.lon), 
                                             drop=True).load();
event_intensity = event_intensity.expand_dims(dim='intensity', axis=3)
events_contour = event.fillna(0)

plt.rc('font', size=12); plt.figure(figsize=(15,6));
for i in enumerate(range(1,len(event.time))):
    ax = plt.subplot(2,3,i[1],projection=ccrs.PlateCarree())
    event_intensity.isel(time=i[0], intensity=0).plot(transform=ccrs.PlateCarree(), vmin=-2, vmax=2,  
                                                      cmap='RdBu_r', extend='both', add_colorbar=True, add_labels=False)
    plt.title(event.isel(time=i[0]).time.values.astype('datetime64[D]'))
    ax.coastlines(resolution='110m', color='black', linewidth=1) 
    ax.add_feature(cfeature.LAND, facecolor='w');
    events_contour.isel(time=i[0]).plot.contour(levels=[34,35], transform=ccrs.PlateCarree(), colors='b', linewidths=4, add_colorbar=False, add_labels=False)

Explore event metrics

from skimage.measure import regionprops_table 

# Use regionprops to calculate it's area and intensity
props = regionprops_table(np.asarray(event.astype('int')),
                          intensity_image=np.asarray(event_intensity.astype('int')),
                          properties=('label', 'area', 'mean_intensity', 'max_intensity'))
props

{'label': array([35]),
 'area': array([31087]),
 'mean_intensity-0': array([0.42535465]),
 'max_intensity-0': array([2])}

Other Examples

%%HTML 
<video controls autoplay width="100%" height="100%" stype="float">
  <source src="./video/GoM_movie.mp4" type="video/mp4">
</video>

%%HTML 
<video controls autoplay width="100%" height="100%" stype="float">
  <source src="./video/blob_movie.mp4" type="video/mp4">
</video>

Global History and Accounting of Marine Heatwaves

What new insights does Ocetrac reveal about marine heatwaves?

• Two different types of marine heatwave patterns emerge: localized events that grow in place and globally connected events that are linked through the tropics


• Marine heatwaves with global connectivity tend to be longer lasting and more intense compared to localized events


• The tropics provide an important conduit for extremely large-scale and persistent marine heatwaves


• Marine heatwaves split and merge – connected in space and time

Why use Ocetrac?

• Overcomes the limitations of characterizing the complex spatial connectivity and temporal behavior of marine heatwaves as they evolve


• Enables the characterization of new spatial patterns and behaviors of some of the most dangerous marine heatwaves of the 21st Century


• Provides a new public catalog of observed marine heatwave metrics and trajectories


• New set of tools to probe other extreme phenomena

Get Involved!

• Try it out

conda install -c conda-forge ocetrac

• File issues and make pull requests on GitHub https://github.com/ocetrac/ocetrac
• E-mail me directly: scannell@ldeo.columbia.edu