SST anomaly / intensity statistics

[j-atkins]

Hello,

First off, thanks for this great tool and comprehensive documentation, I've enjoyed trying it out already! 🙂

I have a question about the event statistics returned by tracker.run(). I see there a lots of nice stats associated with area, centroid, duration etc. but it would be great to also have some on average, cumulative and peak intensities for each identified event/ID. Is that within scope for the project?

It's possible I've missed a good way of doing so in the existing implementation (using tracked ID_field to extract from the original SST fields [anomalies]?) so please do let me know if that's the case and/or the best way to get these kind of stats!

Thanks very much!

Jamie

[wienkers]

Hi Jamie, thanks for the kind words !

The tracker module operates solely on binary fields, so it never sees the originating anomaly data. This was a deliberate design/encapsulation choice that keeps the detection & tracking algorithms more focused & efficient... but it does mean intensity statistics need to be computed as a post-processing step.

These sort of more broadly useful summary intensity stats is something I'd like to add in a future release, likely as a separate pipeline step after tracking. This would still allow the flexibility of studying compound extremes, etc, using whatever anomaly data one wants to feed in.

Until then, you can compute these by combining the tracked ID_field with your anomaly data.
Here's a relatively efficient vectorised approach that computes intensity metrics for all events in parallel, with the result indexed by the ID coordinate:

extremes_ds = xr.open_zarr("extremes_binary.zarr", chunks={})
events_ds = xr.open_zarr("extreme_events.zarr", chunks={})

anomaly = extremes_ds.dat_anomaly
ID_field = events_ds.ID_field
IDs = xr.DataArray(events_ds.ID).chunk(chunks={"ID": 200})  # N.B.: The intermediate arrays would get too large without this ID chunk

# Boolean mask with dimensions (time, lat, lon, ID)
event_mask = ID_field == IDs

# Resulting DataArrays are (time, ID)
mean_intensity = anomaly.where(event_mask).mean(dim=["lat", "lon"])
peak_intensity = anomaly.where(event_mask).max(dim=["lat", "lon"])

# Resulting DataArray is (lat, lon, ID)
cumulative_intensity = anomaly.where(event_mask).sum(dim=["time"])   # We might need to be more clever here, depending on how big your data is...

# Mask where the event isn't present
mean_intensity = mean_intensity.where(events_ds.presence).compute()
peak_intensity = peak_intensity.where(events_ds.presence).compute()

From here you can get overall event-level summaries.

Note that the means here are not area-weighted. If you need (proper) area-weighted mean intensities, you can incorporate grid cell areas into the calculation using apply_ufunc in a similar pattern. Happy to point you in the right direction for that if needed !

Hope that helps !
Cheers,
Aaron

[j-atkins]

Great, thanks very much Aaron for the thorough explanation and post-processing script, super useful!