Long-Term Trends

In previous chapters, we reduced the time dimension using aggregations and rolling averages. Now, we will extract the actual rate of change. We will move beyond visual inspection and apply mathematical modeling to our data cubes to quantify exactly how fast our environment is changing.

Applied Scenario We will continue using the NOAA Extended Reconstructed Sea Surface Temperature dataset. Our goal is to calculate the linear rate of ocean warming over the last century, interpret the mathematical baseline of our models, and ultimately apply harmonic regression to see how the Earth’s seasonal cycle influences our long-term trend estimates.

1. The Mathematics of Change¶

In the previous chapter, we observed that the global sea surface temperature (SST) increase followed an approximately linear trajectory. Visualizing your data is always a crucial first step before committing to a specific statistical trend model.

While we will focus on linear trends, there are several more complex models used in environmental data science depending on the behavior of the data:

• Polynomial Regression: Fits curves (such as quadratic or cubic equations) to data where the rate of change is actively accelerating or decelerating over time.

• Harmonic Regression: Uses sine and cosine functions to model highly cyclical patterns, well suited for raw time series where the seasonal cycle has not been removed.

• Non-parametric Smoothing (e.g., LOESS): Fits localized curves to data to visualize underlying patterns without forcing a single global mathematical equation.

Before we scale up to entire maps, we need to understand the core idea: a trend is the result of fitting a straight line through data over time. When we talk about long-term climate trends, we typically refer to time series of 30 years or more. In this example, we use linear regression as a simple and robust method for quantifying these trends. We want to fit a straight line through our noisy temporal data.

The equation for a straight line is modeled as:

• $y$ (value) is our observed variable, such as temperature.

• $x$ (time) is our time index.

• $b$ (intercept) is the value at the beginning of the time series.

• $m$ (slope) represents how fast the variable changes over time.

Visual Intuition¶

In climate science, the slope is the most critical component. It transforms a complex time series into a single, highly interpretable number describing environmental change.

• Positive slope: An increasing trend (e.g., warming). A steep slope means rapid change.

• Negative slope: A decreasing trend (e.g., cooling).

• Zero slope: A flat line, meaning the variable is stable over time with little to no long-term change.

For an alternative standalone version of the explorer, follow this link.

2. Point Based Trends First¶

Before we run a calculation across the entire globe, it is always best practice to test our statistical model on a single point. This ensures we understand the inputs and outputs before scaling up.

We will select a single coordinate in the North Atlantic and use xarray’s built-in .polyfit() method. The deg=1 argument specifies that we want a first-degree polynomial, which is a straight line.

The output is a new dataset containing polyfit_coefficients. This array has a new dimension called degree. The value at degree=1 is our slope, and the value at degree=0 is our intercept.

Visualizing the point trend¶

To verify that our model worked, we can use xarray.polyval() to evaluate those coefficients and generate the actual trend line, which we can then plot over our raw observations.

3. Pixel Wise Polyfit¶

Now we scale up from a single point to all pixels. Because xarray is designed for multidimensional arrays, scaling up to the entire planet requires no manual loop structures. Instead, we use vectorized operations across the time dimension.

We apply the exact same function to our master dataset, and xarray broadcasts the linear regression across every pixel simultaneously.

To truly appreciate what just happened, let us take an immediate visual peek at the resulting slope array. Even before we convert the mathematical units into something readable, we can map the raw values using a diverging colormap (RdBu_r) and robust=True to automatically handle outliers.

You have just computed tens of thousands of individual linear regression models in a single line of code. This is the power of labeled multidimensional arrays and vectorized operations.

4. Interpreting the Slope¶

When calculating a trend, the slope inherently has units of Variable per Time Step (e.g., °C per time step). But what exactly is a single time step in our dataset?

We must discuss units carefully. When you feed datetime64 coordinate objects into xarray.polyfit(), the underlying algorithm automatically converts the dates into nanoseconds to perform the continuous math. Therefore, the raw slope returned by our function is currently measured in degrees Celsius per nanosecond.

A temperature change per nanosecond is infinitely small and scientifically meaningless to read. We need to scale this slope into a standard, interpretable metric: degrees Celsius per decade.

Now we have a highly interpretable variable. For example, if a pixel has a value of 0.2, it means that specific patch of ocean has been warming at a rate of 0.2°C every ten years.

Visualizing the scaled trend¶

Let’s plot our newly scaled data. Because our values now cross zero (featuring both warming and cooling), we must use a diverging colormap (like RdBu_r). Crucially, we must also pass center=0 to ensure that 0°C/decade is perfectly mapped to white, allowing the red and blue colors to accurately represent warming and cooling.

5. Interpreting the Intercept¶

While the slope tells us how fast temperatures are changing, the intercept ($b$) tells us our starting point. But what exactly is the “start” of our time series mathematically?

When xarray converts datetime64 objects to nanoseconds for the regression, it uses the standard Unix epoch. This means that mathematical time zero ($x=0$) is exactly January 1, 1970. Therefore, the intercept returned by our model represents the modeled sea surface temperature on that specific date.

It is important to distinguish between the start of your specific dataset and the start of the underlying time measurement. By pure chance, the time series we are analyzing might also begin on January 1, 1970. However, even if your dataset started in 1854 or 2000, the mathematical intercept calculated by xarray would still represent the modeled value on January 1, 1970, simply because of how the Unix epoch anchors the timeline.

Let us extract and visualize the intercept (where degree=0). Since these are absolute temperatures rather than diverging rates of change, we will use a sequential colormap like viridis.

Why does this look like a normal temperature map?¶

You might notice that this intercept map closely resembles a standard global distribution of actual sea surface temperatures, with warm equators and cold poles. This is not a coincidence; it highlights a critical concept in climate data.

The geographic variation in temperature across the Earth is enormous, often spanning more than 30°C from the equator to the freezing poles. In contrast, the long term climate trend we extracted earlier is relatively small (measured in fractions of a degree per decade). Because the climate change signal is so small compared to the massive spatial temperature gradient, the mathematical intercept mirrors the physical state of the ocean.

Concept Check: The Unix Epoch Anchor¶

You use .polyfit() on a dataset of monthly temperatures extending from 2000 to 2020. What does the resulting intercept ($b$) represent mathematically?

A) The modeled temperature on January 1, 2000 (the start of your dataset). B) The modeled temperature on January 1, 1970 (the Unix epoch). C) The global average temperature over the entire 20-year period.

6. Modeling the Seasonal Cycle¶

In our previous examples, we relied on a simple linear regression to calculate long term trends. However, we lumped all the data together and ignored the inherent seasonal variations.

What if we want to model the long-term trend while explicitly including these seasonal variations? To capture both a steady trend and a repeating seasonal cycle simultaneously, we use harmonic regression. This model combines a standard straight line with sine and cosine waves.

The Mathematics of the Harmonic Model¶

The harmonic model is fitted to every pixel by estimating four distinct parameters. The equation takes the following form:

• $\beta_0$ (Intercept): The baseline value at the start of the time series.

• $\beta_1$ (Linear Trend): The long term rate of change (e.g., global warming).

• $\beta_2$ and $\beta_3$ (Harmonic Terms): These coefficients control the amplitude and phase of the seasonal cycle.

• $T$ (Period): The length of one full cycle. Because our data follows an annual cycle, $T = 1$ year.

By fitting this model, we compress decades of complex spectral-temporal dynamics into just four interpretable numbers per pixel.

For an alternative standalone version of the explorer, follow this link.

Preparing the Time Coordinate¶

Unlike .polyfit(), which automatically converts dates to nanoseconds, xarray provides a more flexible .curvefit() method for custom models. For .curvefit() to work smoothly and for our period ($T=1$) to make mathematical sense, we need to convert our datetime coordinate into fractional years (e.g., July 1st, 1970 becomes 1970.5).

Point Based Harmonic Fit¶

Just as we did with the linear regression, we will start by fitting this complex model to a single grid cell in the North Atlantic. First, we define our mathematical model as a standard Python function.

Visualizing the Fitted Values¶

To understand how well our model fits the raw data, we can plug our estimated $\beta$ coefficients back into the harmonic_model function to generate the fitted curve.

Scaling to the Global Map¶

Just like .polyfit(), .curvefit() can be broadcast across the entire multidimensional array by calling it directly on ds.sst.

Note: Pixel-wise curve fitting is highly processing intensive. It requires the computer to run an optimization algorithm on tens of thousands of individual time series. If you run this globally, be prepared for it to take significantly longer than the simple linear regression!

Once processed, visualizing the coefficients reveals drastically different patterns. Plotting $\beta_2$ or $\beta_3$ will reveal geographic maps of seasonal amplitude, highlighting regions with extreme summer to winter temperature swings. However, our primary interest is the long-term warming trend, which is captured by $\beta_1$.

Visualizing the Harmonic Slope¶

Before we plot the map, we need to consider our units. This is a great opportunity to see how custom models differ from built-in functions.

Because we manually defined our time coordinate (frac_year) in years, the $\beta_1$ slope is already calculated in degrees Celsius per year. Unlike the nanosecond conversion required for .polyfit(), we only need to multiply by 10 to reach our climate standard of degrees per decade.

7. Comparing the Models¶

We have now calculated the global warming trend using two different methods on the raw dataset:

1. Simple Linear Regression (.polyfit()): Fits a straight line through the raw, seasonally noisy data.

2. Harmonic Regression (.curvefit()): Fits a straight line while simultaneously modeling the seasonal sine and cosine waves.

Which one is better? Let us compute the mathematical difference between the two slopes to see if accounting for seasonality changed our trend estimates.

Interpreting the Difference¶

Looking at the map, you will immediately notice that the difference between the two models is not just random noise. Instead, there is a striking, hemisphere-wide pattern: the Northern Hemisphere is predominantly purple (where the linear model estimated a higher trend), while the Southern Hemisphere is predominantly orange (where the harmonic model estimated a higher trend).

Why do the hemispheres show opposite biases? The answer lies in the seasons.

The Northern and Southern hemispheres experience summer and winter exactly six months out of phase. If a dataset starts or ends during a specific time of year (for example, ending in December when the North is cold and the South is warm), a simple linear regression gets “tricked” by these final data points. It interprets the normal seasonal high or low at the end of the time series as part of the overall trajectory, slightly skewing the long-term slope up or down.

Because the harmonic model explicitly calculates the sine and cosine waves of the seasonal cycle, it is immune to this start- and end-date bias. It mathematically separates the seasonal phase from the long-term trend. This hemisphere-scale map perfectly illustrates why harmonic regression is considered the more rigorous approach when analyzing raw, non-anomalous environmental data.

8. Exercise: Seasonal Trend Analysis¶

Task: Combine the temporal grouping skills you learned in the previous chapter with your new trend-fitting skills. Your goal is to map the long-term warming trend specifically for the month of January.

1. Filter: Select only the January data from the original ds.sst dataset. (Hint: Use .sel() combined with the .dt accessor for logical indexing, or use .groupby("time.month")).

2. Model: Fit a degree 1 polynomial to this January subset.

3. Convert: Extract the slope and convert it from nanoseconds to degrees Celsius per decade.

4. Map: Plot the resulting slope map using a diverging colormap explicitly centered on zero.

5. Interpret: Write two short sentences comparing your January map to the annual average map we created earlier.

Starter code:

Sample solution

# 1. Filter for January (Month 1)
january_data = ds.sst.sel(time=ds.sst.time.dt.month == 1)

# 2. Model the trend
jan_trend = january_data.polyfit(dim="time", deg=1)
jan_slope_raw = jan_trend.polyfit_coefficients.sel(degree=1)

# 3. Convert units (nanoseconds to decades)
ns_per_year = 1e9 * 60 * 60 * 24 * 365.25
jan_slope_decade = (jan_slope_raw * ns_per_year) * 10

# 4. Map the results
plt.figure(figsize=(12, 6))
jan_slope_decade.plot(
    cmap="RdBu_r", 
    center=0, 
    vmin=-0.3, 
    vmax=0.3, 
    cbar_kwargs={"label": "January Warming Rate (°C / Decade)"}
)
plt.title("January Sea Surface Temperature Trend")
plt.show()

Interpretation:

While the January map largely mirrors the broad global warming pattern, it reveals intensified regional cooling in the Southern Ocean and the eastern Pacific. This highlights how specific seasonal extremes, such as austral summer conditions, can drive localized temperature changes that are otherwise muted in a year-round average.

9. Summary¶

By applying statistical models across multidimensional arrays, we can transition from observing temporal patterns to precisely measuring the rate of environmental change.

• Pixel-Wise Modeling: A trend map represents a statistical model computed at every single pixel simultaneously. Using xarray, this is achieved through vectorized operations rather than manual loops.

• Unit Conversion: When using .polyfit() with standard datetime coordinates, the slope is calculated in units per nanosecond and must be manually scaled to a meaningful metric (like degrees per decade).

• Interpreting the Intercept: The mathematical starting point ($x=0$) for xarray datetime regressions is the standard Unix epoch. Consequently, the intercept always represents the modeled baseline value on January 1, 1970.

• Harmonic Regression: While linear regression works well on anomalies, fitting raw data with a harmonic model (.curvefit()) captures both the seasonal cycle and the long-term trend. This prevents seasonal start- and end-date biases from skewing your overall trend estimate.

• Diverging Colormaps: When visualizing rates of change that span both positive (warming) and negative (cooling) values, always use a diverging colormap explicitly centered on zero (center=0).