Space-for-Time Method

🌍 Proof of Concept: Space-for-Time Analysis

Estimate land cover change impacts using spatial variability

Author: Daniel E. Pabon-Moreno

Overview

Land use changes have significant consequences on the Earth system's energy budget. For example, deforestation of tropical rainforests profoundly impacts carbon and water cycles at local, regional, and global scales.

The space-for-time method uses local vegetation contrasts to disentangle the effects of land use changes on biophysical variables like Land Surface Temperature (LST).

What You'll Learn

The theoretical basis of the space-for-time method
How to generate synthetic test data
How to run and interpret space4time_proc results
How to filter results using R² and co-occurrence thresholds

Setup

using Pkg
Pkg.add(url="https://github.com/JuliaStats/GLM.jl", rev="f4047d4930957bc5317fd4d0b73f197383a4ee4a")
Pkg.add(url = "https://github.com/dpabon/YAXArraysToolbox.jl")
using SkipNan
using YAXArraysToolbox
using YAXArrays
using Zarr
using CairoMakie
using GeoMakie
using Random
using NeutralLandscapes
using TiledViews
using DimensionalData
using Statistics

   Resolving package versions...
     Project No packages added to or removed from `~/work/YAXArraysToolbox.jl/YAXArraysToolbox.jl/docs/Project.toml`
    Manifest No packages added to or removed from `~/work/YAXArraysToolbox.jl/YAXArraysToolbox.jl/docs/Manifest.toml`
    Updating git-repo `https://github.com/dpabon/YAXArraysToolbox.jl`
   Resolving package versions...
     Project No packages added to or removed from `~/work/YAXArraysToolbox.jl/YAXArraysToolbox.jl/docs/Project.toml`
    Manifest No packages added to or removed from `~/work/YAXArraysToolbox.jl/YAXArraysToolbox.jl/docs/Manifest.toml`

# Set temporary directory for YAXArrays
YAXArrays.YAXdir("/tmp/YAXA_tmp")

# Activate CairoMakie for plotting
CairoMakie.activate!()

The Space-for-Time Method

Conceptual Framework

The space-for-time method is based on a key assumption: spatial variability in land cover within a local area can serve as a proxy for temporal changes in land cover.

🔑 Key Insight

If we observe how Land Surface Temperature varies across different land cover types within a small moving window, we can estimate the temperature change that would occur if the land cover changed over time.

Algorithm Steps

The method follows these steps:

Define a moving window over the study area
Within each window, analyze the relationship between land cover composition and the target variable
Use regression to estimate the contribution of each land cover type
Calculate transition effects when converting between land cover types

flowchart TD
    A[Land Cover Map] --> B[Moving Window]
    C[Biophysical Variable] --> B
    B --> D[Local Regression]
    D --> E[Land Cover Coefficients]
    E --> F[Transition Effects]

Generating Synthetic Data

To demonstrate the method, we'll create synthetic land cover and temperature data with known properties.

Step 1: Create a Land Cover Map

We'll generate a 1000 × 1000 pixel land cover map (1 pixel = 1 meter, so 1 km × 1 km total):

edge = 1000
size_tile = (edge, edge)

(1000, 1000)

Generate a spatially autocorrelated pattern using midpoint displacement:

Random.seed!(232323)
spatial_auto = 0.9
midpoint_sim = rand(MidpointDisplacement(spatial_auto), size_tile)
heatmap(midpoint_sim)

Classify into 3 land cover classes:

n_classes = 3
classes_dist = NeutralLandscapes.classify(midpoint_sim, ones(n_classes))

1000×1000 Matrix{Float64}:
 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0  …  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0  …  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 ⋮                        ⋮              ⋱            ⋮                   
 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0  …  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0
 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0     1.0  1.0  1.0  1.0  1.0  1.0  1.0

heatmap(classes_dist)

Step 2: Assign LST Values per Class

We define a constant Land Surface Temperature (LST) for each class:

# LST values per class (°C)
class_1_lst = 20.0  # e.g., Forest (coolest)
class_2_lst = 22.0  # e.g., Grassland
class_3_lst = 23.8  # e.g., Cropland (warmest)

all_lst = (class_1_lst, class_2_lst, class_3_lst)

# Create LST map
lst = fill(NaN, size_tile)

for i in eachindex(all_lst)
    lst[findall(==(i), classes_dist)] .= all_lst[i]
end

Spatial Resampling

Land cover maps typically have coarser resolution than 1 meter. We'll aggregate to a new spatial resolution and estimate the frequency of each class per pixel.

Resample Land Cover Classes

size_pixel_new = 20  # New pixel size: 20×20 = 400 original pixels

# Use TiledView to create views of each new pixel
a = TiledView(classes_dist, (size_pixel_new, size_pixel_new), (0, 0); keep_center = false, pad_value = NaN32)

# Array for resampled class frequencies
new_res_array_classes = fill(0.0, (size_pixel_new, size_pixel_new, n_classes))

for i in 1:size_pixel_new
    for j in 1:size_pixel_new
        for c in 1:n_classes
            new_res_array_classes[i, j, c] = count(==(c), a[:, :, i, j]) / (size_pixel_new^2)
        end
    end
end

Visualize the class frequencies:

Class 1:

fig = Figure()
ax = Axis(fig[1, 1]; xlabel = "x", ylabel = "y", title = "Class 1 Frequency")
temp = heatmap!(new_res_array_classes[:, :, 1], colormap = Reverse(:bamako))
Colorbar(fig[1, 2], temp, label = "Occurrence")
fig

Class 2:

fig = Figure()
ax = Axis(fig[1, 1]; xlabel = "x", ylabel = "y", title = "Class 2 Frequency")
temp = heatmap!(new_res_array_classes[:, :, 2], colormap = Reverse(:bamako))
Colorbar(fig[1, 2], temp, label = "Occurrence")
fig

Class 3:

fig = Figure()
ax = Axis(fig[1, 1]; xlabel = "x", ylabel = "y", title = "Class 3 Frequency")
temp = heatmap!(new_res_array_classes[:, :, 3], colormap = Reverse(:bamako))
Colorbar(fig[1, 2], temp, label = "Occurrence")
fig

Resample LST

a = TiledView(lst, (size_pixel_new, size_pixel_new), (0, 0); keep_center = false)

new_res_array_lst = fill(NaN, (size_pixel_new, size_pixel_new))

for i in 1:size_pixel_new
    for j in 1:size_pixel_new
        new_res_array_lst[i, j] = mean(a[:, :, i, j])
    end
end

fig = Figure()
ax = Axis(fig[1, 1]; xlabel = "x", ylabel = "y", title = "Mean LST")
temp = heatmap!(new_res_array_lst[:, :], colormap = :lajolla)
Colorbar(fig[1, 2], temp, label = "LST (°C)")
fig

Organizing Data as YAXArrays

Land Cover Cube

axlist = (
    lon(1:size(new_res_array_classes, 1)),
    lat(1:size(new_res_array_classes, 2)),
    Variables(["class$i" for i in 1:n_classes])
)

lcc_cube = YAXArray(axlist, new_res_array_classes)
lcc_cube

┌ 20×20×3 YAXArray{Float64, 3} ┐
├──────────────────────────────┴──────────────────────────────────── dims ┐
  ↓ lon Sampled{Int64} 1:20 ForwardOrdered Regular Points,
  → lat Sampled{Int64} 1:20 ForwardOrdered Regular Points,
  ↗ Variables Categorical{String} ["class1", …, "class3"] ForwardOrdered
├─────────────────────────────────────────────────────── loaded in memory ┤
  data size: 9.38 KB
└─────────────────────────────────────────────────────────────────────────┘

LST Cube

axlist_lst = (
    lon(1:size(new_res_array_lst, 1)),
    lat(1:size(new_res_array_lst, 2))
)

lst_cube = YAXArray(axlist_lst, new_res_array_lst)
lst_cube

┌ 20×20 YAXArray{Float64, 2} ┐
├────────────────────────────┴──────────────────────── dims ┐
  ↓ lon Sampled{Int64} 1:20 ForwardOrdered Regular Points,
  → lat Sampled{Int64} 1:20 ForwardOrdered Regular Points
├───────────────────────────────────────── loaded in memory ┤
  data size: 3.12 KB
└───────────────────────────────────────────────────────────┘

Running Space4Time Analysis

Now we can use the space4time_proc function:

results = space4time_proc(
    lst_cube,
    lcc_cube;
    time_axis_name = nothing,
    classes_var_name = :Variables,
    classes_vec = ["class1", "class2", "class3"],
    winsize = 5,
    minpxl = 0,
    minDiffPxls = 0,
    max_value = 1,
    showprog = true
)

YAXArray Dataset
Shared Axes: 
  (↓ lon Sampled{Int64} 1:20 ForwardOrdered Regular Points,
  → lat Sampled{Int64} 1:20 ForwardOrdered Regular Points)

Variables with additional axes:
  Additional Axes: 
  (↓ classes Categorical{String} ["class1", …, "class3"] ForwardOrdered,
  → values_of_Z_for_pure_classes Categorical{String} ["estimated", "estimated_error"] ForwardOrdered)
  Variables: 
  metrics_for_classes

  Additional Axes: 
  (↓ summary_stat Categorical{String} ["rsquared_adjusted", …, "predicted"] Unordered)
  Variables: 
  summary_mov_window

  Additional Axes: 
  (↓ transitions Categorical{String} ["class1 to class2", …, "class2 to class3"] ForwardOrdered,
  → differences Categorical{String} ["delta", …, "co_occurrence"] Unordered)
  Variables: 
  metrics_for_transitions

results

YAXArray Dataset
Shared Axes: 
  (↓ lon Sampled{Int64} 1:20 ForwardOrdered Regular Points,
  → lat Sampled{Int64} 1:20 ForwardOrdered Regular Points)

Variables with additional axes:
  Additional Axes: 
  (↓ classes Categorical{String} ["class1", …, "class3"] ForwardOrdered,
  → values_of_Z_for_pure_classes Categorical{String} ["estimated", "estimated_error"] ForwardOrdered)
  Variables: 
  metrics_for_classes

  Additional Axes: 
  (↓ summary_stat Categorical{String} ["rsquared_adjusted", …, "predicted"] Unordered)
  Variables: 
  summary_mov_window

  Additional Axes: 
  (↓ transitions Categorical{String} ["class1 to class2", …, "class2 to class3"] ForwardOrdered,
  → differences Categorical{String} ["delta", …, "co_occurrence"] Unordered)
  Variables: 
  metrics_for_transitions

Interpreting Results

The space4time analysis returns several outputs:

Output	Description
`coefficients`	Estimated effect of each land cover class
`R²`	Goodness of fit for each moving window
`transition effects`	Expected change when converting between classes

Expected Results

Given our synthetic data with known LST values:

Transition	Expected ΔT (°C)	Interpretation
Class 1 → Class 2	+2.0	Forest → Grassland warming
Class 1 → Class 3	+3.8	Forest → Cropland warming
Class 2 → Class 3	+1.8	Grassland → Cropland warming

Filtering Results

Using R² and Co-occurrence Thresholds

metrics_transitions_cube = results.metrics_for_transitions

┌ 3×3×20×20 YAXArray{Union{Missing, Float64}, 4} ┐
├────────────────────────────────────────────────┴─────────────────────── dims ┐
  ↓ transitions Categorical{String} ["class1 to class2", …, "class2 to class3"] ForwardOrdered,
  → differences Categorical{String} ["delta", …, "co_occurrence"] Unordered,
  ↗ lon Sampled{Int64} 1:20 ForwardOrdered Regular Points,
  ⬔ lat Sampled{Int64} 1:20 ForwardOrdered Regular Points
├──────────────────────────────────────────────────────────── loaded in memory ┤
  data size: 28.12 KB
└──────────────────────────────────────────────────────────────────────────────┘

metrics_transitions_cube[differences = At("delta")].data

3×20×20 view(::Array{Union{Missing, Float64}, 4}, :, 1, :, :) with eltype Union{Missing, Float64}:
[:, :, 1] =
 NaN  NaN    NaN    NaN    NaN    NaN    …  NaN    NaN    NaN    NaN    NaN
 NaN  NaN    NaN    NaN    NaN    NaN       NaN    NaN    NaN    NaN    NaN
 NaN    1.8    1.8    1.8    1.8    1.8       1.8    1.8    1.8    1.8    1.8

[:, :, 2] =
 NaN  NaN    NaN    NaN    NaN    NaN    …  NaN    NaN    NaN    NaN    NaN
 NaN  NaN    NaN    NaN    NaN    NaN       NaN    NaN    NaN    NaN    NaN
 NaN    1.8    1.8    1.8    1.8    1.8       1.8    1.8    1.8    1.8    1.8

[:, :, 3] =
 NaN  NaN    NaN    NaN    NaN    NaN    …  NaN    NaN    NaN    NaN    NaN
 NaN  NaN    NaN    NaN    NaN    NaN       NaN    NaN    NaN    NaN    NaN
 NaN    1.8    1.8    1.8    1.8    1.8       1.8    1.8    1.8    1.8    1.8

;;; … 

[:, :, 18] =
   2.0    2.0    2.0    2.0    2.0  …  2.0  2.0  2.0  2.0  NaN    NaN
 NaN    NaN    NaN    NaN    NaN       3.8  3.8  3.8  3.8  NaN    NaN
 NaN    NaN    NaN    NaN    NaN       1.8  1.8  1.8  1.8    1.8    1.8

[:, :, 19] =
   2.0    2.0    2.0    2.0    2.0  …  2.0  2.0  2.0  2.0  NaN    NaN
 NaN    NaN    NaN    NaN    NaN       3.8  3.8  3.8  3.8  NaN    NaN
 NaN    NaN    NaN    NaN    NaN       1.8  1.8  1.8  1.8    1.8    1.8

[:, :, 20] =
   2.0    2.0    2.0    2.0    2.0  …    2.0  2.0  2.0  2.0  NaN    NaN
 NaN    NaN    NaN    NaN    NaN       NaN    3.8  3.8  3.8  NaN    NaN
 NaN    NaN    NaN    NaN    NaN       NaN    1.8  1.8  1.8    1.8    1.8

Apply quality filters:

masking_without_delta = masking_proc(
    results.metrics_for_transitions;
    cube_rsquared = results.summary_mov_window[summary_stat = At("rsquared_adjusted")],
    rsquared_thr = 0.2,
    cube_co_occurrence = results.metrics_for_transitions[Differences = At("co_occurrence")],
    co_occurence_thr = 0.5,
    cube_delta = nothing,
    time_dim = nothing,
    showprog = true
)

┌ 3×3×20×20 YAXArray{Union{Missing, Float64}, 4} ┐
├────────────────────────────────────────────────┴─────────────────────── dims ┐
  ↓ transitions Categorical{String} ["class1 to class2", …, "class2 to class3"] ForwardOrdered,
  → differences Categorical{String} ["delta", …, "co_occurrence"] Unordered,
  ↗ lon Sampled{Int64} 1:20 ForwardOrdered Regular Points,
  ⬔ lat Sampled{Int64} 1:20 ForwardOrdered Regular Points
├──────────────────────────────────────────────────────────── loaded in memory ┤
  data size: 28.12 KB
└──────────────────────────────────────────────────────────────────────────────┘

Plotting Results

Compare the estimated deltas with the original (known) values:

delta_1_org = abs(class_1_lst - class_2_lst)
delta_2_org = abs(class_1_lst - class_3_lst)
delta_3_org = abs(class_2_lst - class_3_lst)

1.8000000000000007

vec_delta_orig = Array{Float64}(reshape(
    results.metrics_for_transitions[differences = At("delta")].data,
    (size_pixel_new^2 * n_classes)
))

vec_delta_plot = vec_delta_orig[findall(!isnan, vec_delta_orig)]

vec_index = repeat(1:n_classes, outer = (size_pixel_new^2 * n_classes))
vec_index = vec_index[findall(!isnan, vec_delta_orig)]

501-element Vector{Int64}:
 3
 3
 3
 3
 3
 3
 3
 3
 3
 3
 ⋮
 3
 1
 2
 3
 1
 2
 3
 3
 3

transitions = lookup(masking_without_delta, :transitions)
fig = Figure()
ax = Axis(fig[1,1], xticks = (1:length(transitions), transitions))
temp = scatter!(ax, vec_index, vec_delta_plot)
temp2 = scatter!(ax, [i + 0.1 for i in 1:n_classes],
    [delta_1_org, delta_2_org, delta_3_org],
    marker = :diamond, color = :red, markersize = 15)
Legend(fig[1, 2],
    [temp, temp2],
    ["Space4time Results", "Original Delta"])
fig

Validation

The space4time method successfully recovers the known temperature differences between land cover classes, demonstrating the validity of the approach.

Summary

✅ Key Takeaways

The space-for-time method estimates land cover change impacts without long time series
Spatial variability within moving windows serves as a proxy for temporal change
Quality filtering using R² and co-occurrence improves result reliability
The method can quantify biophysical effects of land use transitions

References

Duveiller, Gregory, Josh Hooker, and Alessandro Cescatti. "A Dataset Mapping the Potential Biophysical Effects of Vegetation Cover Change." Scientific Data 5, no. 1 (2018): 1. https://doi.org/10.1038/sdata.2018.14
Li, Yan, Maosheng Zhao, Safa Motesharrei, Qiaozhen Mu, Eugenia Kalnay, and Shuangcheng Li. "Local Cooling and Warming Effects of Forests Based on Satellite Observations." Nature Communications 6, no. 1 (2015): 6603. https://doi.org/10.1038/ncomms7603

Next Steps

📖 Return to Basic Operations for core function usage
📚 See the API Reference for complete function documentation