Spatial Detrend

Note

Go to the end to download the full example code.

Spatial Detrend#

Before proceeding with all the steps, first import some necessary libraries and packages

import xarray as xr
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import signal
import easyclimate as ecl

The significance of removing spatial trends in meteorological data is profound, as it allows scientists to isolate and analyze the specific atmospheric phenomena they are interested in without the obfuscating influence of large-scale, systematic patterns. In essence, it is a fundamental pre-processing step that transforms raw data into a form primed for deeper scientific inquiry.

To illustrate its significance, we will first generate sample data. coords

time = pd.date_range("2000-01-01", periods=50, freq="YS")
lat = np.linspace(-90, 90, 20)
lon = np.linspace(-180, 180, 30)

lat_da = xr.DataArray(lat, dims="lat", coords={"lat": lat})
lon_da = xr.DataArray(lon, dims="lon", coords={"lon": lon})

# deterministic components
base = 15.0
i = xr.DataArray(np.arange(len(time)), dims="time", coords={"time": time})
slope = 0.1
time_trend = slope * i

lat_pattern = -0.5 * abs(lat_da) / 90
lon_pattern = 0.01 * (lon_da / 180)

# raw noise
np.random.seed(42)
noise_raw = xr.DataArray(
    np.random.normal(0, 0.5, (len(time), len(lat), len(lon))),
    dims=["time", "lat", "lon"],
    coords={"time": time, "lat": lat, "lon": lon},
)

# make noise have NO linear component along time (per grid point)
noise_detr = noise_raw.reduce(signal.detrend, dim="time")  # removes (a*i+b)

# build data with a perfectly linear trend + detrended noise
data_trend = base + lat_pattern + lon_pattern + time_trend + noise_detr

easyclimate.calc_detrend_spatial can help us remove temporal trends from spatial data.

detr = ecl.calc_detrend_spatial(data_trend, "time")
detr

<xarray.DataArray (lat: 20, lon: 30, time: 50)> Size: 240kB
array([[[ 0.27663096,  0.40429774,  0.08594522, ..., -0.9404822 ,
         -0.6943026 ,  0.52427975],
        [ 0.11030256, -0.28688591, -0.04574397, ...,  0.50976759,
         -0.60785443, -0.15052855],
        [ 0.24386127,  0.35408779, -0.02029858, ..., -1.0707691 ,
          0.07443379, -0.31605801],
        ...,
        [-0.01605039, -0.25980946, -0.26343359, ..., -0.02953023,
          0.19407818,  1.13637225],
        [-0.1860738 ,  0.14253493,  0.59703897, ..., -0.44836873,
         -0.20887582, -0.52824839],
        [-0.15975767,  0.53452703,  1.09403483, ..., -0.44890051,
         -0.21480025,  0.62982466]],

       [[-0.08685433, -0.63800667, -0.07629357, ...,  0.50659556,
          0.56221022, -0.90430023],
        [ 0.96165415,  0.79988727, -0.65019233, ..., -0.3462656 ,
         -0.55894048,  0.91229701],
        [ 0.00329822, -0.07070987, -0.0375997 , ..., -0.27751378,
          0.30549162, -0.64661266],
...
        [-0.55034765,  0.23778265, -0.12998224, ...,  0.23227481,
         -0.36162164, -0.49836971],
        [-0.93294527,  0.30137388,  0.06714513, ...,  0.57604787,
          0.14870593,  0.41452806],
        [-0.43942524,  0.1709705 ,  0.02793663, ..., -0.00505752,
         -0.50785852,  0.45859735]],

       [[ 0.23813895,  0.60129156,  0.23652598, ...,  0.69639076,
         -0.06031024,  0.21204329],
        [ 0.83102362,  0.13248527,  0.83849968, ...,  0.12382099,
         -0.58467096,  0.1617484 ],
        [ 0.17967965,  0.06744452, -0.29685271, ..., -0.2692266 ,
         -0.17270445,  0.44583756],
        ...,
        [-0.28831745,  0.46964097, -0.21228535, ..., -0.49678879,
         -0.51662595, -0.4253716 ],
        [-0.18999564,  0.14234134,  0.56662112, ...,  0.33723332,
         -0.52528493, -1.00059652],
        [ 0.06775525,  0.57146426, -0.08966109, ...,  0.6495117 ,
         -1.20656517,  0.26882577]]], shape=(20, 30, 50))
Coordinates:
  * lat      (lat) float64 160B -90.0 -80.53 -71.05 -61.58 ... 71.05 80.53 90.0
  * lon      (lon) float64 240B -180.0 -167.6 -155.2 ... 155.2 167.6 180.0
  * time     (time) datetime64[us] 400B 2000-01-01 2001-01-01 ... 2049-01-01


Next, we consider the correct validation: detr(data_trend) should equal noise_detr.

err = np.max(np.abs(detr - noise_detr))
print("max abs error (should be ~1e-12~1e-10):", float(err))

max abs error (should be ~1e-12~1e-10): 1.4932499681208355e-14

Here we also compared the difference before and after removing the trend at a specific point, which can be said to align with our expectations.

data_trend.isel(lat = 10, lon = 20).plot(label = "Trend data")
detr.isel(lat = 10, lon = 20).plot(label = "Detrend data")
plt.legend()
lat = 4.737, lon = 68.28
<matplotlib.legend.Legend object at 0x7f268547fc10>

To better handle large spatial datasets, easyclimate also provides easyclimate.calc_detrend_spatial_fast implementations for significantly faster solutions. These methods include: "scipy_reduce","scipy","numpy","rust","rust_chunked","rust_flexible". Their results demonstrate considerable robustness.

detr = ecl.calc_detrend_spatial_fast(data_trend, "time", method = "scipy_reduce")
# ✅ Correct validation: detr(data_trend) should equal noise_detr
err = np.max(np.abs(detr - noise_detr))
print("[scipy_reduce] max abs error (should be ~1e-12~1e-10):", float(err))

detr = ecl.calc_detrend_spatial_fast(data_trend, "time", method = "scipy")
# ✅ Correct validation: detr(data_trend) should equal noise_detr
err = np.max(np.abs(detr - noise_detr))
print("[scipy] max abs error (should be ~1e-12~1e-10):", float(err))

detr = ecl.calc_detrend_spatial_fast(data_trend, "time", method = "rust")
# ✅ Correct validation: detr(data_trend) should equal noise_detr
err = np.max(np.abs(detr - noise_detr))
print("[rust] max abs error (should be ~1e-12~1e-10):", float(err))

detr = ecl.calc_detrend_spatial_fast(data_trend, "time", method = "rust_chunked")
# ✅ Correct validation: detr(data_trend) should equal noise_detr
err = np.max(np.abs(detr - noise_detr))
print("[rust_chunked] max abs error (should be ~1e-12~1e-10):", float(err))

detr = ecl.calc_detrend_spatial_fast(data_trend, "time", method = "rust_flexible")
# ✅ Correct validation: detr(data_trend) should equal noise_detr
err = np.max(np.abs(detr - noise_detr))
print("[rust_flexible] max abs error (should be ~1e-12~1e-10):", float(err))

[scipy_reduce] max abs error (should be ~1e-12~1e-10): 1.4932499681208355e-14
[scipy] max abs error (should be ~1e-12~1e-10): 1.4932499681208355e-14
[rust] max abs error (should be ~1e-12~1e-10): 1.8596235662471372e-14
[rust_chunked] max abs error (should be ~1e-12~1e-10): 1.8596235662471372e-14
[rust_flexible] max abs error (should be ~1e-12~1e-10): 1.8596235662471372e-14

Total running time of the script: (0 minutes 0.527 seconds)