Spatial Data Analysis in Google Earth Engine

A Geospatial Analyst’s Workflow, Steps 1–6: Acquisition, Preprocessing, Exploratory Analysis, Association, Regression, and Validation

1 Why Spatial Analysis Needs Its Own Workflow

If you already work with tabular data, you know the arc: clean the data, visualize it, test associations, model it. Spatial analysis follows the same arc, but at every stage you are carrying an extra piece of information that tabular data doesn’t have: where. That “where” is what makes spatial analysis harder than it first looks, and it is also what makes it powerful.

Two ideas sit underneath everything that follows, and it’s worth stating them plainly before touching any code:

  • Tobler’s First Law of Geography: “Everything is related to everything else, but near things are more related than distant things.” This is why we can’t always treat spatial observations as independent — a core assumption behind ordinary regression — and it’s why spatial statistics exists as its own discipline.
  • The Modifiable Areal Unit Problem (MAUP): the boundaries you choose to aggregate data into (e.g., wards vs. sub-counties) can change your results, even when the underlying data hasn’t changed. Keep this in mind whenever you pick a zonal unit later in this workflow.

This document walks through all six stages of a spatial workflow: Steps 1–3 run in Google Earth Engine (GEE), a cloud platform that gives you API access to petabytes of satellite imagery and geospatial datasets without needing to download anything locally (code shown in the GEE JavaScript API, the Code Editor at code.earthengine.google.com). Steps 4–6 run in R, because formal spatial statistics — bivariate spatial association, spatial regression, and rigorous validation — are not something GEE implements; GEE is the imagery and aggregation engine, R is the statistics engine, and the export at the end of Step 2 is the handoff point between them.

2 Key Terms Before We Start

Term Definition
AOI (Area of Interest) The geographic boundary — point, line, or polygon — that defines the spatial extent of your analysis. Everything you compute is clipped or filtered to this.
CRS (Coordinate Reference System) The mathematical framework that ties coordinates (e.g., latitude/longitude) to actual locations on Earth’s surface. Two datasets in different CRSs will not overlay correctly until reprojected to a common one.
Projection The specific method of flattening the Earth’s curved surface onto a 2D plane (e.g., UTM, Albers Equal Area). Choice of projection affects area, distance, and shape accuracy — not just appearance.
Spatial resolution (pixel size / scale) The ground area represented by a single pixel, e.g., Sentinel-2 at 10m means each pixel covers a 10m × 10m square on the ground.
Band A single layer of a satellite image capturing reflectance in one part of the electromagnetic spectrum (e.g., Red, Near-Infrared). Multiple bands stacked together form the full image.
ImageCollection GEE’s data structure for a stack of images over time (e.g., every Sentinel-2 scene over your AOI for a year).
FeatureCollection GEE’s data structure for a set of vector geometries with attributes (e.g., administrative boundaries, sample points).
Cloud masking The process of flagging and excluding cloud- or shadow-contaminated pixels from analysis, using quality-assurance bands provided with the imagery.
Compositing / Mosaicking Combining multiple images (e.g., all cloud-free scenes in a month) into a single representative image, usually via a per-pixel statistic like median or max.
Spectral index A formula combining two or more bands to isolate a specific physical signal — e.g., NDVI isolates vegetation greenness from raw reflectance.
Resampling Recalculating pixel values when changing a raster’s resolution or alignment, using a method such as nearest-neighbor (for categorical data) or bilinear/cubic (for continuous data).
Reducer A GEE function that collapses many pixel values into a summary statistic — e.g., ee.Reducer.mean() — either across time (compositing) or across space (zonal statistics).
Zonal statistics Summarizing pixel values within defined polygons (e.g., mean NDVI per ward). The spatial equivalent of “group by” in tabular analysis.
Mixed pixel / edge effect Distortion introduced when a pixel spans two zones or land-cover types, so its value is not purely representative of either.
Spatial autocorrelation The degree to which nearby locations have more similar values than distant ones (see Tobler’s Law). It must be checked before applying ordinary (non-spatial) regression, because strong autocorrelation violates the independence assumption behind standard statistical tests.
Moran’s I A statistic measuring global spatial autocorrelation, ranging roughly from −1 (perfect dispersion) through 0 (random/no pattern) to +1 (perfect clustering). It answers: “overall, are similar values located near each other more than chance would predict?”
Spatial weights matrix A definition of which locations count as “neighbors” of each other (e.g., shared boundary, distance threshold, k-nearest-neighbors), and how strongly each neighbor relationship is weighted. Every spatial autocorrelation statistic depends on this choice.
LISA (Local Indicators of Spatial Association) A local, per-zone version of Moran’s I — instead of one global statistic, each zone gets its own value showing whether it sits inside a local cluster or is a spatial outlier.
Hotspot analysis (Getis-Ord Gi*) A statistic identifying statistically significant spatial clusters of high values (“hot spots”) or low values (“cold spots”), as distinct from clusters of similar values in either direction (which is what Moran’s I/LISA detect).
Choropleth map A map that shades polygons (zones) according to the value of a variable — the standard way to visualize zonal statistics.
Join count statistic A test for spatial autocorrelation in categorical (not continuous) spatial data — e.g., does a binary land-use class tend to cluster next to itself more than chance would predict? The categorical-data counterpart to Moran’s I.
Lee’s L statistic A bivariate statistic measuring spatial association between two different variables while explicitly accounting for spatial autocorrelation in each — the spatial analog of a correlation coefficient, but correcting for the fact that nearby observations are not independent.
OLS (Ordinary Least Squares) Standard, non-spatial linear regression. Valid only when residuals are independent — a condition spatial autocorrelation routinely violates, which is why OLS residuals must be checked with Moran’s I before trusting OLS results on spatial data.
Spatial Lag Model (SAR lag) A regression model that adds a spatially lagged version of the outcome variable as a predictor, for cases where the outcome itself diffuses across space (e.g., disease spread between neighboring areas).
Spatial Error Model (SEM) A regression model that accounts for spatial autocorrelation in the error term rather than the outcome — appropriate when clustering arises from unmeasured spatially structured variables (e.g., unobserved soil quality affecting a regional outcome).
Lagrange Multiplier (LM) test A diagnostic test applied to OLS residuals that tells you whether a Spatial Lag Model or a Spatial Error Model (or neither) is the statistically better-justified choice.
GWR (Geographically Weighted Regression) A regression technique that lets coefficients vary smoothly across space, rather than forcing one global coefficient on the whole study area — it directly models spatial non-stationarity (the idea that a relationship, e.g. rainfall → NDVI, may be stronger in some places than others).
Bandwidth (in GWR) The spatial extent (distance or number of neighbors) used to calculate each local regression in GWR. Too small → overfit, noisy local estimates. Too large → GWR collapses toward a global OLS result. Usually chosen by cross-validation or corrected AIC (AICc).
AICc (corrected Akaike Information Criterion) A model-comparison statistic, penalized for sample size, commonly used to select GWR bandwidth and to compare competing spatial models. Lower is better.
Spatial cross-validation A validation strategy that splits data into spatial blocks (rather than random rows) before train/test splitting, to prevent nearby, spatially autocorrelated points from leaking information between the training and test sets and inflating apparent accuracy.
RMSE (Root Mean Square Error) A standard accuracy metric for continuous predictions — the average magnitude of prediction error, in the same units as the outcome variable.
Confusion matrix For categorical/classification outputs (e.g., land-cover classes), a table cross-tabulating predicted vs. actual classes, from which accuracy metrics (overall accuracy, kappa, producer’s/user’s accuracy) are derived.

3 Step 1: Data Acquisition & Inspection

Purpose of this stage: get the right data, for the right place, at the right time, and confirm — before you do anything else — that it is trustworthy and correctly aligned. This is the spatial equivalent of “cleaning” in tabular analysis, but because spatial data is layered (multiple bands, multiple time steps, multiple sources with different native grids), inspection has to happen at several levels, not just one.

3.1 1.1 Define the Area of Interest (AOI)

What this step does: establishes the geographic boundary that scopes every later computation. In GEE this is either a ee.Geometry (a shape you define directly) or an ee.FeatureCollection (a set of polygons pulled from an existing boundary dataset, such as administrative units).

Analyst’s note: This is the single most consequential decision in the whole workflow, because every later filter, mosaic, and statistic inherits whatever is wrong with the AOI. Two practical rules:

  • Prefer an authoritative boundary dataset (FAO GAUL, GADM, national statistics office shapefiles) over hand-drawn polygons, unless your question specifically requires a custom extent.
  • Check the AOI’s CRS. Administrative boundary datasets are usually in geographic coordinates (EPSG:4326, i.e., latitude/longitude in decimal degrees) — fine for filtering and display, but not for accurate area or distance calculations, which need a projected CRS (see §1.6).
// Option A: Use a known administrative boundary dataset
var aoi = ee.FeatureCollection("FAO/GAUL/2015/level2")
  .filter(ee.Filter.eq('ADM1_NAME', 'Kakamega'));

// Option B: Define a custom polygon directly
var aoiGeom = ee.Geometry.Polygon([[
  [34.70, 0.10], [34.95, 0.10], [34.95, 0.35], [34.70, 0.35], [34.70, 0.10]
]]);

Map.centerObject(aoi, 9);
Map.addLayer(aoi, {color: 'red'}, 'AOI boundary');

3.2 1.2 Search and load the relevant imagery or feature collection

What this step does: loads the raw satellite imagery or dataset you’ll work with, filtered to your AOI and date range.

Analyst’s note — choosing a sensor is a trade-off, not a default:

Sensor / dataset Native resolution Revisit frequency Typical use
Sentinel-2 10–20m ~5 days Fine-scale vegetation, land cover, small-area change
Landsat 8/9 30m 16 days Long time-series (back to 1980s for Landsat family), moderate-scale change
MODIS 250m–1km Daily Regional/continental trend analysis, not property-scale work
CHIRPS ~5.5km Daily Rainfall estimation over large areas

The resolution and revisit frequency you need should be driven by the scale of the phenomenon you’re studying — don’t default to the highest-resolution sensor available if your question is regional.

var s2 = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
  .filterBounds(aoi)
  .filterDate('2024-01-01', '2024-12-31');

print('Number of images found:', s2.size());
print('Metadata of first image:', s2.first());

3.3 1.3 Filter by cloud cover and quality metadata

What this step does: removes whole scenes that are too cloud-contaminated to be useful, using scene-level metadata (as opposed to §1.4, which masks individual pixels within scenes you keep).

var s2Filtered = s2.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20));
print('Images remaining after cloud filter:', s2Filtered.size());

Analyst’s note: Always print .size() immediately before and after every filter you apply. An empty collection is the most common silent failure in GEE — a .mosaic() or .median() call on zero images does not throw an error; it just quietly produces a blank layer, and you won’t notice until a later step behaves strangely. Treat .size() checks as a mandatory habit, not an optional debugging step.

3.4 1.4 Mask individual cloud/shadow pixels

What this step does: even scenes that pass the scene-level cloud filter in §1.3 usually still have some cloudy or cirrus-contaminated pixels. Cloud masking flags those specific pixels as “no data” using a quality-assurance (QA) band delivered with the image, rather than discarding the whole scene.

function maskS2clouds(image) {
  var qa = image.select('QA60');
  var cloudBitMask = 1 << 10;   // bit 10 flags cloud
  var cirrusBitMask = 1 << 11;  // bit 11 flags cirrus
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
    .and(qa.bitwiseAnd(cirrusBitMask).eq(0));
  return image.updateMask(mask).divide(10000)
    .copyProperties(image, ['system:time_start']);
}

var s2Masked = s2Filtered.map(maskS2clouds);

Analyst’s note: This is the direct spatial equivalent of handling missing or erroneous values in a spreadsheet — except “missing” here is defined per-pixel, per-scene, and must be resolved before you combine multiple scenes into one composite, or the cloud contamination will bleed into your final surface.

3.5 1.5 Composite (mosaic) multiple scenes into one analysis-ready image

What this step does: collapses a stack of images over a time period into a single representative image, using a per-pixel reducer (see Key Terms).

var composite = s2Masked.median().clip(aoi);

Map.addLayer(composite, {
  bands: ['B4', 'B3', 'B2'],  // Red, Green, Blue → natural color
  min: 0, max: 0.3
}, 'True-color composite');

Analyst’s note — know what your reducer choice implies:

  • Median is the robust default: it suppresses residual cloud/shadow noise and outlier pixels better than a mean would. Good for a general “typical conditions over the period” surface.
  • Mean is sensitive to outliers (a missed cloud pixel will skew it).
  • Max (e.g., for NDVI) captures peak greenness — useful for phenology studies, but it is not a “typical conditions” surface.

A median composite over a year will blur out short-lived events (a flood, a fire scar). If your question is event-specific, use a tightly date-filtered single scene instead of compositing.

3.6 1.6 Verify projection, resolution, and CRS consistency

What this step does: confirms the coordinate system and pixel size of your working image before you combine it with any other dataset.

print('Projection:', composite.select('B4').projection());
print('Native scale (meters):', composite.select('B4').projection().nominalScale());

Analyst’s note: This is the step most often skipped, and the one most likely to quietly invalidate a multi-source analysis later. Three concrete risks:

  1. Mismatched CRS. If your AOI boundary is in EPSG:4326 and your imagery is delivered in a UTM zone, GEE will reproject on the fly during operations like reduceRegion — but silently, and not always in the way you’d choose deliberately. Decide your working CRS explicitly.
  2. Mismatched resolution. Combining 10m Sentinel-2 data with 5.5km CHIRPS rainfall data requires an explicit decision about which resolution to standardize on — this is handled in Step 2.
  3. Area distortion in geographic CRS. EPSG:4326 is not equal-area — a degree of longitude covers less ground distance near the poles than at the equator. Any area calculation (e.g., “hectares of forest lost”) must be done in a projected, equal-area CRS, not directly in EPSG:4326.

4 Step 2: Preprocessing & Derived Variables

Purpose of this stage: transform the clean, aligned imagery from Step 1 into the specific variables your analysis actually needs, and make sure any second data source you bring in (e.g., rainfall, population, soil type) is placed on a common, defensible spatial grid before you compare it to anything.

4.1 2.1 Compute spectral indices (spatial feature engineering)

What this step does: combines raw reflectance bands mathematically to isolate a specific physical or biological signal. This is the spatial analog of engineered features in tabular modeling (e.g., a BMI computed from height and weight).

// NDVI (Normalized Difference Vegetation Index) — vegetation greenness/health
// Formula: (NIR - Red) / (NIR + Red)
var ndvi = composite.normalizedDifference(['B8', 'B4']).rename('NDVI');

// NDWI (Normalized Difference Water Index) — surface water presence
// Formula: (Green - NIR) / (Green + NIR)
var ndwi = composite.normalizedDifference(['B3', 'B8']).rename('NDWI');

// EVI (Enhanced Vegetation Index) — vegetation signal with reduced
// atmospheric and soil-background noise, better in high-biomass areas
var evi = composite.expression(
  '2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))', {
    'NIR': composite.select('B8'),
    'RED': composite.select('B4'),
    'BLUE': composite.select('B2')
}).rename('EVI');

var indices = composite.addBands([ndvi, ndwi, evi]);

Analyst’s note: NDVI values run from −1 to 1. Roughly: water and bare soil sit near 0 or negative, sparse vegetation sits around 0.2–0.4, and dense healthy vegetation sits above 0.6. Always sanity-check your computed index against these expected ranges before trusting downstream statistics — a badly masked cloud pixel often shows up as an implausible NDVI value.

4.2 2.2 Bring in a second data source and reconcile its resolution

What this step does: loads an additional variable (here, rainfall) and explicitly resamples it to a common working resolution with the imagery from Step 1, so that a pixel-by-pixel or zone-by-zone comparison is arithmetically valid.

var rainfall = ee.ImageCollection("UCSB-CHG/CHIRPS/DAILY")
  .filterBounds(aoi)
  .filterDate('2024-01-01', '2024-12-31')
  .sum()  // total annual rainfall, mm
  .clip(aoi);

// Explicitly resample and reproject to a defined common grid
var rainfallResampled = rainfall.resample('bilinear')
  .reproject({crs: 'EPSG:4326', scale: 30});

Analyst’s note: Resampling CHIRPS from ~5.5km down to 30m does not create new rainfall information — it does not make the estimate more spatially precise than the original data actually was. It only makes the pixel grids align so a comparison with 30m NDVI is computable. Say this explicitly in any write-up: a reader glancing at a smooth 30m rainfall map could easily over-interpret its precision if you don’t flag the resampling.

Choosing a resampling method: - Nearest-neighbor — required for categorical data (e.g., land-cover class codes), because averaging category codes produces meaningless values. - Bilinear / cubic — appropriate for continuous data (rainfall, temperature, NDVI), because interpolating between neighboring values is physically meaningful.

4.3 2.3 Zonal statistics — aggregate pixel values within administrative or sampling units

What this step does: summarizes pixel-level values into per-zone statistics — this is the direct spatial equivalent of a “group by” + aggregate operation in tabular analysis, where the grouping variable is a polygon rather than a categorical column.

var zonalStats = indices.select('NDVI').reduceRegions({
  collection: aoi,       // FeatureCollection of sub-units (e.g., wards)
  reducer: ee.Reducer.mean().combine({
    reducer2: ee.Reducer.stdDev(),
    sharedInputs: true
  }),
  scale: 10,
  tileScale: 4            // raise this if you hit memory/computation errors
});

print('Zonal NDVI statistics (first 5 zones):', zonalStats.limit(5));

Analyst’s note — this is where MAUP (see Introduction) becomes concrete. The zonal mean you get for “NDVI per ward” depends on how those ward boundaries were drawn. If your zones vary hugely in size or shape, consider reporting both the mean and the standard deviation (as above) so a reader can see how much within-zone variability the aggregate is hiding.

4.4 2.4 Handle edge effects and mixed pixels

What this step does: accounts for the distortion introduced when pixels straddle a zone boundary or contain more than one land-cover type.

Analyst’s note — three practical mitigations:

  1. Set an explicit tileScale (4–8) on large reduceRegions calls — this is a computational safeguard, but it also forces you to think about whether your AOI is too large for the resolution you’ve chosen.
  2. ee.Reducer.mean() inside reduceRegions is area-weighted by default, meaning a pixel that is 40% inside a polygon contributes 40% of its value — confirm this is actually the behavior you want for your reducer choice, especially with ee.Reducer.mode() or count-based reducers.
  3. For small polygons relative to pixel size (e.g., 10m pixels over very small wards), consider whether zonal statistics are even reliable — below a certain zone-to-pixel size ratio, mixed-pixel error can dominate the estimate.

4.5 2.5 Export analysis-ready data for downstream statistics

What this step does: hands off the processed, aligned data to a file format usable outside GEE.

Export.table.toDrive({
  collection: zonalStats,
  description: 'NDVI_zonal_stats_2024',
  fileFormat: 'CSV'
});

Export.image.toDrive({
  image: indices,
  description: 'Indices_composite_2024',
  scale: 10,
  region: aoi.geometry(),
  maxPixels: 1e13
});

Analyst’s note: GEE is strong through exactly this point — acquisition, masking, index computation, and zonal aggregation, all at planetary scale without local storage. For formal spatial statistics (Moran’s I test for spatial autocorrelation, Geographically Weighted Regression, spatial lag/ error models), export here and continue in R (spdep, spgwr, mgwr) or Python (pysal), which is where your existing regression and association toolkit resumes.

5 Step 3: Exploratory Spatial Data Analysis (ESDA)

Purpose of this stage: this is the spatial equivalent of plotting histograms, boxplots, and a correlation matrix before modeling tabular data — except here the first and most important “plot” is the map itself, and the first and most important statistical check is whether your data is spatially clustered at all. GEE handles visualization natively; the formal autocorrelation statistics (Moran’s I, Getis-Ord Gi*) require exporting your zonal data (from §2.3/§2.5) and computing them in R or Python, because GEE does not implement spatial-weights-based statistics itself. This split is intentional and worth stating plainly to a reader: GEE is the imagery and aggregation engine; R/Python is the statistics engine.

5.1 3.1 Visualize the zonal variable as a choropleth map

What this step does: renders your zonal statistic (e.g., mean NDVI per ward, from §2.3) as a shaded map — the single fastest way to spot obvious spatial patterns, data errors, and candidate clusters before running any formal test.

// Visualize the zonal NDVI mean computed in Step 2.3
var ndviVis = zonalStats.style({
  color: 'black',
  fillColor: '00000000',
  width: 1
});

// Create a choropleth by painting each zone with its NDVI_mean value
var choropleth = ee.Image().byte().paint({
  featureCollection: zonalStats,
  color: 'NDVI_mean'
});

Map.addLayer(choropleth, {
  min: 0, max: 0.8,
  palette: ['white', 'yellow', 'green', 'darkgreen']
}, 'Mean NDVI by zone');

Analyst’s note: Choose your color palette and class breaks deliberately — a poorly chosen color ramp (e.g., a rainbow palette with no perceptual ordering) can visually exaggerate or hide real patterns. A sequential palette (light → dark, as above) is the correct default for a continuous variable like NDVI; reserve diverging palettes (e.g., blue–white–red) for variables with a meaningful midpoint, such as change from a baseline.

5.2 3.2 Export zonal data and construct a spatial weights matrix

What this step does: moves the zonal table (already exported in §2.5) into R or Python and defines which zones count as “neighbors” of each other — a decision every subsequent spatial statistic depends on.

# R — using sf and spdep
library(sf)
library(spdep)

zones <- st_read("NDVI_zonal_stats_2024.shp")

# Queen contiguity: zones sharing any boundary point are neighbors
nb <- poly2nb(zones, queen = TRUE)
weights <- nb2listw(nb, style = "W")  # row-standardized weights

Analyst’s note: The choice of neighbor definition is not a technicality — it materially changes your results. Queen contiguity (shares an edge or a corner) is a reasonable default for administrative polygons. Distance-based or k-nearest-neighbor weights are better suited to point data (e.g., health-facility locations) where “contiguity” doesn’t apply. State your choice explicitly in any methods section; a reviewer or collaborator should be able to reproduce it.

5.3 3.3 Test for global spatial autocorrelation (Moran’s I)

What this step does: answers the question “is this variable spatially clustered overall, or could its spatial pattern plausibly be random?” — this is the spatial equivalent of checking a correlation coefficient before building a regression model, and the result directly determines whether ordinary (non-spatial) regression is even appropriate later.

moran_test <- moran.test(zones$NDVI_mean, weights)
print(moran_test)

Analyst’s note — how to read the output:

  • Moran’s I near 0 → little to no spatial clustering; values appear roughly randomly distributed across zones.
  • Moran’s I significantly > 0 (check the p-value) → positive spatial autocorrelation: similar values (high-high or low-low) cluster together. This is the typical finding for environmental variables like NDVI.
  • Moran’s I significantly < 0 → negative spatial autocorrelation: a checkerboard-like pattern where neighboring zones differ sharply. This is less common but does occur (e.g., strongly zoned land use).

If Moran’s I is significant, you now have statistical justification for why an ordinary regression later would understate your standard errors — the observations are not spatially independent, so you’ll need a spatial regression approach (spatial lag or spatial error model, or GWR) when you get there.

5.4 3.4 Identify local clusters and outliers (LISA)

What this step does: breaks the single global Moran’s I statistic down to the zone level, so you can see where the clustering is happening and classify each zone into one of four types: High-High (a high value surrounded by high values), Low-Low, High-Low (a high-value outlier among low neighbors), or Low-High.

local_moran <- localmoran(zones$NDVI_mean, weights)
zones$lisa_cluster <- attr(local_moran, "quadr")$mean

# Map the cluster classification
plot(zones["lisa_cluster"], main = "LISA Cluster Map: NDVI")

Analyst’s note: A LISA map is often more useful to a decision-maker than the global Moran’s I alone — “which specific wards are low-vegetation outliers surrounded by high-vegetation neighbors” is an actionable finding, whereas “the region shows significant clustering overall” is not, by itself, something a programme officer can act on.

5.5 3.5 Hotspot analysis (Getis-Ord Gi*)

What this step does: identifies statistically significant clusters specifically of high values (hot spots) or low values (cold spots), which is a distinct question from LISA’s “similar values near each other” — Gi* tells you whether a cluster’s local sum is significantly higher or lower than you’d expect if values were randomly distributed in space.

library(spdep)

gi_star <- localG(zones$NDVI_mean, weights)
zones$gi_star_z <- as.numeric(gi_star)

# Z-scores beyond ±1.96 are significant at the 95% confidence level
zones$hotspot_class <- cut(zones$gi_star_z,
  breaks = c(-Inf, -1.96, 1.96, Inf),
  labels = c("Cold spot", "Not significant", "Hot spot"))

Analyst’s note: Moran’s I/LISA and Getis-Ord Gi* are complementary, not redundant — run both. LISA tells you about similarity-based clustering (including low-low clusters, which Gi* on its own tends to under-emphasize), while Gi* gives a cleaner, more interpretable hot-spot/cold-spot map for a non-technical audience, such as a ministry of health partner deciding where to target an intervention.

6 Step 4: Spatial Association Analysis

Purpose of this stage: this is the spatial equivalent of correlation and chi-square testing in tabular analysis — but a plain Pearson correlation or chi-square test run on spatial data without adjustment is misleading, because spatially autocorrelated observations are not independent, and both tests assume independence. Step 4 shows the corrected, spatially-aware versions of these familiar tools.

We continue directly from the zones object and weights spatial weights matrix built in §3.2, now assuming zones holds two variables of interest — NDVI_mean (from Step 2) and rain_mean (zonal mean of the resampled rainfall from §2.2, extracted the same way as §2.3).

6.1 4.1 Naive (non-spatial) correlation — and why it can mislead

What this step does: computes an ordinary Pearson or Spearman correlation between two zonal variables, as a baseline.

cor.test(zones$NDVI_mean, zones$rain_mean, method = "pearson")

Analyst’s note: Report this number only alongside its spatial correction (§4.2), never alone. If both NDVI and rainfall are positively spatially autocorrelated (which is typical — nearby wards tend to resemble each other in both vegetation and rainfall), a naive correlation test’s p-value will be overstated — it will appear more statistically significant than it really is, because the effective number of independent observations is lower than the raw zone count suggests. This is a direct, practical consequence of Tobler’s Law from the Introduction.

6.2 4.2 Spatially-corrected bivariate association (Lee’s L)

What this step does: computes Lee’s L, a bivariate statistic that measures association between two variables while explicitly incorporating the spatial weights matrix, giving a more honest significance test than §4.1.

library(spdep)

lee_result <- lee.test(zones$NDVI_mean, zones$rain_mean,
                        listw = weights, zero.policy = TRUE)
print(lee_result)

Analyst’s note — how to read Lee’s L alongside the naive correlation:

  • If Lee’s L and the naive Pearson correlation are similar in magnitude and both significant, your finding is robust to spatial structure.
  • If the naive correlation is significant but Lee’s L is not (or is much weaker), the apparent association was likely an artifact of both variables independently clustering in space — not genuine covariation. This distinction matters enormously when the finding will inform a resource-allocation decision.

6.3 4.3 Association between two categorical spatial variables (join counts)

What this step does: when your variables of interest are categorical rather than continuous (e.g., a binary land-degradation flag, or presence/ absence of a disease outbreak per zone), join count statistics test whether zones of the same category cluster next to each other more than chance would predict — the spatial, categorical-data counterpart to a chi-square test of association.

# Example: a binary "high NDVI" vs "low NDVI" classification per zone
zones$ndvi_class <- factor(ifelse(zones$NDVI_mean > median(zones$NDVI_mean),
                                   "High", "Low"))

joincount.test(zones$ndvi_class, weights)

Analyst’s note: A significant join-count result for “High-High” joins tells you that high-NDVI zones are non-randomly adjacent to other high-NDVI zones — formal statistical confirmation of what the choropleth map in §3.1 suggested visually. Pair this with the LISA cluster map from §3.4, which identifies which zones form those clusters.

6.4 4.4 Visualize the bivariate relationship

What this step does: plots the zonal relationship directly, which should always accompany a correlation statistic — a single correlation coefficient can hide a nonlinear relationship or an influential outlier zone.

library(ggplot2)

ggplot(zones, aes(x = rain_mean, y = NDVI_mean)) +
  geom_point(size = 2, alpha = 0.7) +
  geom_smooth(method = "lm", se = TRUE, color = "darkgreen") +
  labs(x = "Mean annual rainfall (mm)", y = "Mean NDVI",
       title = "Zonal NDVI vs. rainfall") +
  theme_minimal()

Analyst’s note: Label or flag any zone that visibly departs from the general trend before moving to regression — an outlier zone at this stage often turns out to be a zonal-statistics artifact (e.g., a very small zone dominated by mixed pixels, §2.4) rather than a genuine deviation worth modeling.

7 Step 5: Spatial Regression Modeling

Purpose of this stage: this is the spatial equivalent of your regression stage — but because Step 3 established that spatial autocorrelation is present (via Moran’s I), an ordinary linear regression (OLS) run on zonal data violates the independence assumption behind its own standard errors and significance tests. Step 5 shows how to diagnose this formally and choose the correct spatial regression model.

7.1 5.1 Fit OLS first — as a diagnostic baseline, not a final model

What this step does: fits an ordinary regression, which you will use only to test its residuals for leftover spatial structure — not to report as your final result.

ols_model <- lm(NDVI_mean ~ rain_mean, data = zones)
summary(ols_model)

Analyst’s note: Do not report OLS coefficients and p-values as your final finding on spatial data without first completing §5.2. This is the single most common statistical error in applied spatial analysis — treating zones as if they were independent sampling units.

7.2 5.2 Test OLS residuals for spatial autocorrelation

What this step does: runs Moran’s I on the OLS residuals specifically (not the raw outcome variable, which you already tested in §3.3) — this tells you whether the regression itself leaves spatial structure unexplained, which is the direct justification for moving to a spatial model.

zones$ols_resid <- residuals(ols_model)
moran.test(zones$ols_resid, weights)

Analyst’s note: If this Moran’s I is not significant, your OLS residuals are behaving like independent noise and OLS may be defensible as your final model — document this check either way, since it is exactly the evidence a reviewer will ask for.

7.3 5.3 Decide between Spatial Lag and Spatial Error models (LM tests)

What this step does: runs Lagrange Multiplier diagnostic tests, which are specifically designed to tell you which spatial model — lag or error — is the better-justified correction, rather than choosing one by assumption.

lm_tests <- lm.LMtests(ols_model, weights, test = "all")
summary(lm_tests)

Analyst’s note — the standard decision rule:

  • If LM-lag is significant and LM-error is not → fit a Spatial Lag Model (§5.4).
  • If LM-error is significant and LM-lag is not → fit a Spatial Error Model (§5.5).
  • If both are significant, check the robust versions of each test (lm.LMtests(..., test = "RLMlag") / "RLMerr") and follow whichever robust test remains significant — this resolves cases where the two diagnostics are “fighting” due to shared variance.

7.4 5.4 Fit a Spatial Lag Model

What this step does: fits a regression where the outcome in each zone is modeled partly as a function of the outcome in neighboring zones — appropriate when the phenomenon itself diffuses across space.

library(spatialreg)

lag_model <- lagsarlm(NDVI_mean ~ rain_mean, data = zones, listw = weights)
summary(lag_model)

Analyst’s note: The key additional parameter here is rho (ρ), the spatial lag coefficient. A significant, positive rho means a zone’s NDVI is partly explained by its neighbors’ NDVI, independent of rainfall — i.e., there is genuine spatial spillover in the outcome itself, not just shared exposure to similar rainfall.

7.5 5.5 Fit a Spatial Error Model

What this step does: fits a regression that absorbs spatial autocorrelation into the error structure — appropriate when clustering is driven by an unmeasured, spatially structured variable rather than diffusion of the outcome itself.

error_model <- errorsarlm(NDVI_mean ~ rain_mean, data = zones, listw = weights)
summary(error_model)

Analyst’s note: The key parameter here is lambda (λ), the spatial error coefficient. A significant lambda means there’s a spatially structured factor you haven’t measured (e.g., soil type, unrecorded land management practices) driving the residual clustering — worth naming explicitly as a limitation, since it signals a variable missing from the model rather than a flaw in the modeling approach itself.

7.6 5.6 Compare model fit (AIC) and finalize

What this step does: compares OLS, the lag model, and the error model on a common footing, so your final model choice is evidence-based rather than a default.

AIC(ols_model, lag_model, error_model)

Analyst’s note: Lower AIC indicates better fit relative to model complexity. In almost all cases where Step 3 found significant spatial autocorrelation, both spatial models will outperform OLS — the AIC comparison is there to choose between lag and error, and to document, in your methods section, exactly why OLS was rejected.

7.7 5.7 Account for local non-stationarity: Geographically Weighted Regression

What this step does: so far, §5.4–5.5 still assume the rainfall → NDVI relationship has one fixed strength across your entire AOI. GWR relaxes that assumption and lets the coefficient vary smoothly from zone to zone, revealing where the relationship is strong or weak — a distinct question from “is there spatial autocorrelation,” which the lag/ error models already addressed.

library(GWmodel)

zones_sp <- as(zones, "Spatial")  # GWmodel requires a Spatial object

# Select an optimal bandwidth by corrected AIC
bw <- bw.gwr(NDVI_mean ~ rain_mean, data = zones_sp,
             approach = "AICc", kernel = "bisquare", adaptive = TRUE)

gwr_model <- gwr.basic(NDVI_mean ~ rain_mean, data = zones_sp,
                        bw = bw, kernel = "bisquare", adaptive = TRUE)
print(gwr_model)

Analyst’s note: Map the local coefficient surface (gwr_model$SDF$rain_mean) as a choropleth, the same way you mapped NDVI in §3.1. A coefficient that is strongly positive in one part of your AOI and near zero in another is a substantive finding in its own right — it tells a partner organization that a rainfall-based intervention logic that works in one sub-region may simply not apply in another. This is the kind of result a single global regression coefficient (§5.4–5.6) cannot show you.

8 Step 6: Validation and Accuracy Assessment

Purpose of this stage: this closes the loop — every model fitted in Step 5 needs to be checked against reality before you trust or report it. In spatial analysis, validation has an extra failure mode beyond ordinary train/test splitting: because nearby observations are correlated, a random split can leak information between training and test sets and make your model look more accurate than it actually is.

8.1 6.1 Re-check residuals of your final spatial model

What this step does: confirms that your chosen spatial model (lag, error, or GWR) has actually resolved the autocorrelation problem diagnosed in §5.2 — the final diagnostic check before trusting the model.

zones$final_resid <- residuals(lag_model)  # or error_model, as chosen
moran.test(zones$final_resid, weights)

Analyst’s note: You want this Moran’s I to be non-significant. If it’s still significant after fitting a spatial model, the chosen model hasn’t fully captured the spatial structure — reconsider your spatial weights definition (§3.2) or whether an additional covariate is needed, rather than reporting the model as final.

8.2 6.2 Spatial cross-validation (not a random split)

What this step does: evaluates predictive accuracy using spatially blocked folds, so that training and test observations are not drawn from the same local neighborhood — a random split would let spatially autocorrelated neighbors “leak” into both sets and artificially inflate accuracy.

library(blockCV)
library(sf)

# Create spatial folds (blocks), sized larger than the observed
# spatial-autocorrelation range from Step 3
sb <- cv_spatial(x = zones, column = "NDVI_mean",
                  k = 5, size = 20000,  # block size in meters — set from
                                        # your AOI's autocorrelation range
                  selection = "random")

# Loop over folds, refitting the model on training blocks and
# predicting on the held-out block each time (sketch)

Analyst’s note: Set the block size from evidence, not a round number — a reasonable starting point is the distance range over which Moran’s I (§3.3) remains significant (a variogram, if you have one, gives this more precisely). A block that’s too small defeats the purpose and still leaks spatial information between folds.

8.3 6.3 Accuracy metrics appropriate to your outcome type

What this step does: quantifies predictive performance using the metric that matches your outcome’s data type.

# Continuous outcome (e.g., predicted vs. observed NDVI)
predicted <- predict(lag_model)
actual <- zones$NDVI_mean

rmse <- sqrt(mean((predicted - actual)^2))
r_squared <- cor(predicted, actual)^2
print(c(RMSE = rmse, R2 = r_squared))

# Categorical outcome (e.g., classified land-cover vs. ground-truth points)
# library(caret)
# confusionMatrix(predicted_class, reference_class)

Analyst’s note: Report RMSE in the outcome’s natural units (e.g., “NDVI units” or “mm of rainfall”) so a non-technical reader can judge whether the error is practically meaningful, not just statistically small. For classification accuracy, always report the kappa statistic alongside overall accuracy — overall accuracy alone can look deceptively high on an imbalanced class distribution (e.g., mostly one land-cover type).

8.4 6.4 Ground-truth validation against independent field data

What this step does: compares model output against data that played no role in fitting the model — ideally independently collected field observations (GPS points, survey records) rather than another remote-sensed product, which would just be comparing one model’s assumptions to another’s.

ground_truth <- st_read("field_validation_points.shp")

validation <- st_join(ground_truth, zones["NDVI_mean"])
cor.test(validation$ndvi_observed, validation$NDVI_mean)

Analyst’s note: This is the step that most distinguishes a rigorous spatial analysis from a plausible-looking map. A model can pass every internal diagnostic in §6.1–6.3 and still misrepresent ground conditions if the underlying imagery, index, or zonal aggregation had a systematic bias — independent field validation is the only check that catches that.

9 Summary Checklist

Step Task What it accomplishes
1.1 Define AOI Sets the geographic scope for everything downstream
1.2 Load imagery/collection Selects data source by resolution/revisit trade-off
1.3 Filter by cloud cover Removes unusable whole scenes
1.4 Mask cloud/shadow pixels Handles “missing data” at the pixel level
1.5 Composite/mosaic Produces one representative analysis surface
1.6 Verify CRS/resolution Prevents silent misalignment before combining data
2.1 Compute spectral indices Feature engineering from raw reflectance bands
2.2 Align a second data source Makes cross-source comparison arithmetically valid
2.3 Zonal statistics Spatial “group by” summarization per unit
2.4 Handle edge/mixed pixels Controls bias from boundary and small-zone effects
2.5 Export Hands off to R/Python for formal spatial statistics
3.1 Choropleth visualization Fastest way to spot patterns and data errors before testing
3.2 Build spatial weights matrix Defines “neighbor” relationships that every later statistic depends on
3.3 Global Moran’s I Tests whether the variable is spatially clustered at all
3.4 Local LISA clusters Locates where clustering/outliers occur, zone by zone
3.5 Hotspot analysis (Gi*) Identifies significant high- or low-value clusters for targeting decisions
4.1 Naive correlation Baseline association, flagged as potentially overstated
4.2 Lee’s L Spatially-corrected bivariate association — the honest version of 4.1
4.3 Join count statistics Tests association between two categorical spatial variables
4.4 Bivariate scatterplot Visual check for nonlinearity or outlier zones before regression
5.1 OLS (diagnostic only) Baseline regression, not a final result on spatial data
5.2 Moran’s I on OLS residuals Confirms whether OLS leaves spatial structure unexplained
5.3 Lagrange Multiplier tests Decides between Spatial Lag and Spatial Error models
5.4 Spatial Lag Model Models outcome diffusion between neighboring zones
5.5 Spatial Error Model Absorbs clustering from unmeasured spatial factors
5.6 AIC comparison Evidence-based final model selection
5.7 GWR Reveals where a relationship is locally strong or weak
6.1 Residual re-check Confirms the final model resolved spatial autocorrelation
6.2 Spatial cross-validation Prevents inflated accuracy from spatially leaking folds
6.3 Accuracy metrics (RMSE/kappa) Quantifies predictive performance appropriately by outcome type
6.4 Ground-truth validation Checks the whole pipeline against independent field data

This document now spans the full workflow, from an empty AOI definition in §1.1 through a field-validated spatial model in §6.4. The thread running through every step is the same one named in the Introduction: spatial observations are not independent, and each stage — aggregation, testing, modeling, and validation — has to account for that explicitly rather than borrowing tabular-data methods unchanged.