Spatial Data Analysis in Google Earth Engine

A Geospatial Analyst’s Workflow — Step 1: Data Acquisition & Inspection | Step 2: Preprocessing & Derived Variables

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 the first two stages of a spatial workflow 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 is shown in the GEE JavaScript API (the Code Editor at code.earthengine.google.com), with Python (geemap / earthengine-api) noted where the workflow differs.

---

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.

---

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 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

What comes next (not covered in this document): exploratory spatial data analysis — testing for spatial autocorrelation with Moran’s I, hotspot/cluster detection (Getis-Ord Gi*) — followed by spatial regression (GWR, spatial lag/error models) and accuracy validation against ground-truth points.