Finding Stable Pixels and Choosing PIAs

Finding Stable Pixels and Choosing PIAs#

This notebook measures how much the reflectance of every pixel has moved across a whole Sentinel-2 archive, and uses that map to choose pseudo-invariant areas (PIAs, also called PIFs): the polygons used in relative radiometric normalization, where a scene is adjusted to match a reference by fitting a per-band regression on surfaces assumed not to have changed.

The measurement is the short part — three lines of code. Most of the notebook is about reading it correctly, because the obvious statistics fail in instructive ways.

Why the bands, and not an index#

Vegetation indices are ratios, so they cancel out multiplicative changes in brightness. A pixel can hold exactly the same NDVI for eight years while its reflectance drifts by a third. That drift is precisely what a PIA must not have, so the analysis runs on the bands themselves, which ndvi2gif exposes as ordinary index names:

index='blue' | 'green' | 'red' | 'nir' | 'swir1' | 'swir2'
index='red_edge1' | 'red_edge2' | 'red_edge3'   # Sentinel-2 only

They come back as surface reflectance in 0–1 on every sensor, so a threshold picked on Sentinel-2 keeps its meaning on Landsat.

What gets measured#

Everything here is a property of the pixel with respect to itself. Nothing is a percentile, a rank or a comparison against the rest of the scene, which means the numbers survive a change of extent, of date, or of sensor:

what it answers

std

how much the reflectance moved, in reflectance units

cv

the same, relative to the pixel’s own level

range

the largest excursion it ever had

drift

how much of that movement went in one direction

count

how many observations back the answer

1. Setup#

import warnings
warnings.filterwarnings('ignore', category=UserWarning, module='geemap.conversion')

import ee
import geemap
from ndvi2gif import NdviSeasonality

# ee.Authenticate()          # only the first time
ee.Initialize(project='your-project-id')
print('✓ Earth Engine initialised')

2. Region of interest#

The example runs over the working extent of the Doñana flood-mapping protocol (SW Spain): roughly 114 × 103 km covering the marshes, the Atlantic coast, the Guadalquivir estuary, several reservoirs and the city of Seville.

Any ROI accepted by NdviSeasonality works: a shapefile, a GeoJSON, a drawn geometry, a DEIMS site id or a Sentinel-2 tile code.

UTM = 'EPSG:32629'
START, END = 2018, 2025

# Doñana protocol extent, in its native projected CRS
roi = ee.Geometry.Rectangle(
    [676440, 4059960, 790380, 4163340], proj=UTM, geodesic=False
)

Map = geemap.Map()
Map.centerObject(roi, zoom=8)
Map.addLayer(roi, {'color': 'red'}, 'Study extent')
Map

3. The dispersion map#

NdviSeasonality composites inside a year: periods splits the year into seasons, months or a single 1 January–31 December window, but never spans several years. To get one raster for the whole archive you reduce the underlying collection directly. The object exposes it as ndvi_col — already filtered to the ROI and cloud-masked — and the index function as d[index]:

def series(index, start=START, end=END, max_cloud_cover=60):
    """Every valid observation of `index` over the full period, as one collection."""
    ns = NdviSeasonality(
        roi=roi, sat='S2', start_year=start, end_year=end, index=index,
        cloud_filter=True, max_cloud_cover=max_cloud_cover, scl_mask=True,
    )
    return ns.ndvi_col.filterDate(f'{start}-01-01', f'{end + 1}-01-01').map(ns.d[index])


col = series('swir1')

std = col.reduce(ee.Reducer.stdDev()).rename('std')
mean = col.mean().rename('mean')
cv = std.divide(mean).rename('cv')
rng = col.max().subtract(col.min()).rename('range')
n_obs = col.reduce(ee.Reducer.count()).rename('n_obs')

That is the whole measurement. std is in reflectance units, so 0.02 means the pixel has held its SWIR1 reflectance to within a couple of percentage points across eight years of acquisitions, in every season.

Note

scl_mask=True applies per-pixel cloud and shadow masking from the Scene Classification Layer. Residual cloud is the easiest way to fake instability, so it is not optional here. max_cloud_cover=60 is deliberately loose in return: the scene-level filter discards whole images, and with per-pixel masking in place a partly cloudy scene still contributes its clear half. Dispersion statistics want observations.

Map = geemap.Map()
Map.centerObject(roi, zoom=8)
Map.addLayer(std, {'min': 0.01, 'max': 0.12,
                   'palette': ['1a9850', 'ffffbf', 'd73027']}, 'std (swir1)')
Map.addLayer(cv, {'min': 0, 'max': 1.5,
                  'palette': ['1a9850', 'ffffbf', 'd73027']}, 'cv (swir1)', False)
Map.addLayer(rng, {'min': 0, 'max': 0.5,
                   'palette': ['1a9850', 'ffffbf', 'd73027']}, 'range (swir1)', False)
Map.addLayer(n_obs, {'min': 100, 'max': 900,
                     'palette': ['440154', '21918c', 'fde725']}, 'valid observations', False)
Map

3.1 Which statistic to trust#

They are not interchangeable, and each fails in its own way:

  • std is the workhorse. It averages over hundreds of observations, so no single bad one moves it much.

  • range is max - min, so it is decided by one observation. A cloud edge the SCL mask did not catch blows up the range of an otherwise perfect pixel. Read it as “something happened to this pixel once”, not as a measure of stability.

  • cv is std divided by the pixel’s own level, which is what makes brightness comparable — and what makes it explode on water, where the denominator goes to zero.

  • count is the quality layer for all of them. A standard deviation from three images means nothing. It also reveals the Sentinel-2 tiling: pixels in the overlap between two MGRS tiles are seen about twice as often.

Use std and the drift of section 5 as the two numbers you decide with; keep cv and range as diagnostics for when something does not add up. Two interpretable numbers beat one composite index that fires without telling you why.

3.2 Choosing the time window#

Worth a minute of measurement rather than a guess. key='count' with periods=1 gives the valid observations per pixel of each year, and over this extent it says:

year

scenes

distinct dates

valid obs/pixel (median)

2015

28

9

—

2016

70

28

—

2017

308

89

75

2018

510

142

91

2019

525

145

105

2020–2024

~520

~145

~100

The archive reaches back to 2015, so the limit is not the catalogue. It is Sentinel-2B, launched in March 2017 and operational around the middle of that year: before it, revisit was 10 days instead of 5. That is why 2017 has 89 distinct dates against the ~145 of a normal year, and why they concentrate in the second half.

So 2017 is usable — 75 observations per pixel is plenty for an annual mean — but it is unevenly sampled within the year, which biases the within-year term of section 4 towards whichever seasons were better covered. A reasonable compromise is to include 2017 in the between-year and drift analysis, where each year contributes a single mean, and leave it out of the within-year term. 2016, with 28 dates, is not usable either way.

3.3 Export it and go look at it#

Over an extent this size a direct download will time out, so the raster goes out as an Earth Engine task. Reflectance-like bands are scaled by 10000 and stored as int16, which keeps the file a quarter of the size of float32 and loses nothing you can see.

dispersion = ee.Image.cat([
    std.multiply(10000),
    cv.multiply(10000),
    rng.multiply(10000),
    mean.multiply(10000),
    n_obs,
]).toInt16().rename(['std', 'cv', 'range', 'mean', 'n_obs']).clip(roi)

task = ee.batch.Export.image.toDrive(
    image=dispersion, description='swir1_dispersion',
    folder='ndvi2gif', region=roi, scale=10, crs=UTM, maxPixels=1e13,
)
# task.start()

In QGIS, load std with a green-to-red stretch over roughly 200–1200 (remember the ×10000) and the quiet surfaces come straight out: runways, quarry floors, dune fields, industrial roofs, open water.

4. Two levels: within the year and between years#

The single std pools two very different kinds of movement, and separating them costs almost nothing:

how it is computed

what it captures

within-year

periods=1, key='variance', averaged over years

scatter inside each year: phenology, tides, illumination, residual cloud

between-year

periods=1, key='mean', then reduced across the collection

movement between years: land-cover change, a reservoir filling and emptying

periods=1 is what makes this work: each composite spans the whole year, so the statistic covers every valid observation of that year, winter and summer together. Seasonality is not noise to be averaged away here — a surface that changes with the season is not invariant, and a PIA has to hold regardless of when the scene being normalized was acquired. Keeping the whole year inside one period puts the seasonal cycle into the within-year term, where it belongs.

def annual(index, key):
    """One composite per year, each covering the full year."""
    ns = NdviSeasonality(
        roi=roi, sat='S2', periods=1, start_year=START, end_year=END,
        index=index, key=key, cloud_filter=True, max_cloud_cover=60, scl_mask=True,
    )
    return ns.get_year_composite().select('p1')


band = 'swir1'

# Variances add, so the split is done on variances and square-rooted at the end
within_var = annual(band, 'variance').mean()              # scatter inside a year
annual_level = annual(band, 'mean')                       # one value per year
between_var = annual_level.reduce(ee.Reducer.variance())  # movement across years

total = within_var.add(between_var).sqrt().rename('sigma')
frac_between = between_var.divide(within_var.add(between_var)).rename('frac_between')

# The same idea in more readable units: brightest year minus darkest
between_range = annual_level.max().subtract(annual_level.min())

Map = geemap.Map()
Map.centerObject(roi, zoom=8)
Map.addLayer(within_var.sqrt(), {'min': 0.01, 'max': 0.12,
                                 'palette': ['1a9850', 'ffffbf', 'd73027']}, 'within-year std')
Map.addLayer(between_var.sqrt(), {'min': 0.005, 'max': 0.06,
                                  'palette': ['1a9850', 'ffffbf', 'd73027']}, 'between-year std')
Map.addLayer(frac_between, {'min': 0, 'max': 0.5,
                            'palette': ['ffffff', '4575b4', '313695']},
             'fraction of variance that is interannual', False)
Map

The two terms reconstruct the total dispersion of section 3 almost exactly: over a test window inside this extent, total and the direct std correlate at 0.9998, and their 5th, 50th and 95th percentiles agree to four decimal places (0.0181 / 0.0239 / 0.0679).

So the split buys no precision — it buys diagnosis. Green in the first layer means “this pixel held its reflectance all year”. Blue in the last one, wherever the others are not green, says the instability is movement between years rather than a wobble inside them. Reservoirs and burnt areas light up; irrigated crops do not, because their variance is seasonal and repeats.

5. Drift: is the pixel going somewhere?#

Dispersion says how much a pixel moves. It does not say whether the movement goes anywhere, and only one of those is fatal:

  • A pixel that oscillates and returns has a high std, but it averages out. Its value in a scene from 2019 and one from 2025 has the same expectation.

  • A pixel that drifts has a value that depends on the date. Since the reference scene has one fixed date, that pixel biases the regression by an amount that grows with how far the scene being normalized sits from it. This one does not average away.

SpatialTrendAnalyzer answers the second question, and periods=1, key='mean' already produced the annual series it needs.

from ndvi2gif.timeseries import SpatialTrendAnalyzer

ns_level = NdviSeasonality(
    roi=roi, sat='S2', periods=1, start_year=START, end_year=END,
    index=band, key='mean', cloud_filter=True, max_cloud_cover=60, scl_mask=True,
)

# 'mann_kendall' returns Sen's slope and intercept, the accumulated magnitude,
# and Kendall's tau: the rank correlation with time, scale-free, from -1 to +1
trend = SpatialTrendAnalyzer(ns_level).calculate_pixel_trends(
    method='mann_kendall', min_observations=5,
)

drift = trend.select('magnitude').abs().rename('drift')   # reflectance over the series
tau = trend.select('tau')                                 # how monotonic the movement is

Map = geemap.Map()
Map.centerObject(roi, zoom=8)
Map.addLayer(trend.select('slope'), {'min': -0.01, 'max': 0.01,
                                     'palette': ['313695', 'ffffff', 'a50026']},
             "Sen's slope (reflectance / year)")
Map.addLayer(drift, {'min': 0, 'max': 0.06,
                     'palette': ['1a9850', 'ffffbf', 'd73027']}, 'accumulated drift')
Map.addLayer(tau, {'min': -1, 'max': 1,
                   'palette': ['313695', 'ffffff', 'a50026']}, "Kendall's tau", False)
Map

Warning

ee.Reducer.kendallsCorrelation() also advertises a p-value band, and ndvi2gif does not return it: it comes back fully masked at every series length tested, including trends scipy.stats.kendalltau scores below 1e-20, and a fully masked band silently masks whatever it is combined with. Judge significance from tau and the number of years.

6. Reading the map: what the reference classes say#

The Doñana protocol has a hand-drawn set of pseudo-invariant areas, built and revised over years of operational use. Measuring the statistics inside them shows what the map means in practice — and why no single number settles the question.

class

brightness

σ (red)

σ (nir)

σ (swir1)

CV worst band

% var. interannual

drift (swir1)

Arena (dunes)

0.457

0.044

0.047

0.059

0.12

0.06

0.006

Fosfoyesos (gypsum stacks)

0.474

0.082

0.081

0.082

0.21

0.21

0.050

Urbano-2

0.263

0.045

0.050

0.057

0.23

0.04

0.004

Aeropuertos (runways)

0.269

0.046

0.046

0.053

0.24

0.09

0.004

Urbano-1

0.260

0.047

0.052

0.061

0.24

0.04

0.006

Minería (quarries)

0.257

0.057

0.058

0.067

0.29

0.13

0.034

Marisma mareal (tidal marsh)

0.178

0.039

0.040

0.040

0.39

0.07

0.003

Pastizales (turf)

0.295

0.047

0.080

0.052

0.58

0.06

0.015

Pinar (pine forest)

0.137

0.032

0.036

0.044

0.65

0.17

0.033

Mar (open sea)

0.018

0.021

0.021

0.022

1.39

0.02

0.000

Embalses (reservoirs)

0.034

0.039

0.046

0.044

1.84

0.08

0.001

Sentinel-2 2018–2025, averaged inside each class of a reference layer in its state of 2026-08-24 (118 polygons, 11 classes), drawn by hand and independently of this method.

Four things in that table are worth internalising before drawing any polygon.

A std threshold favours dark surfaces. Dark pixels have little signal, so they have little variance — not because they are invariant but because there is nothing there to vary. At σ ≤ 0.05, reservoirs (0.039–0.046) pass and dunes (0.044–0.059) do not, which is exactly backwards.

A cv threshold favours bright ones. By CV, Mar and Embalses are the two worst classes in the list, yet open sea is one of the most reliable normalization targets there is. Its CV is an artefact of dividing by 0.018.

Neither of them sees drift. Embalses and Fosfoyesos sit at comparable σ, but the reservoirs have essentially zero drift — they filled, emptied and filled again, and came back — while the gypsum stacks climb 0.008 reflectance per year, 0.050 across the series. Those are the phosphogypsum ponds at Huelva, under active capping and revegetation. Same dispersion, opposite verdicts: a reservoir may be salvageable by picking better polygons, a surface under construction is not.

Look at the bands separately. Pastizales here is irrigated turf — golf greens and football pitches — and it is the only class whose bands disagree: 0.080 in the near infrared against 0.047 in red and 0.052 in SWIR1. Red saturates against chlorophyll and sits still, while the near infrared follows leaf area, which moves with every mowing and with the autumn overseeding. If your normalization is fitted per band, that class is usable in red and SWIR and not in NIR.

Note

Pinar drifting 0.033 with tau 0.39 deserves a look rather than a shrug. The 2017 Las Peñuelas fire burned pine forest at the edge of this extent, and post-fire recovery is exactly the monotonic SWIR1 signal this test is built to find.

7. Drawing the PIAs#

The regression that normalizes a scene needs targets spread across the brightness range: a set clustered at one end constrains slope and offset very poorly. Nothing in the map guarantees that spread, and — as section 6 shows — no single statistic can, because every threshold leans one way or the other.

The class structure is what provides it. Dark sea, dark-to-mid pine and marsh, mid urban fabric and runways, bright dunes: the classes already stratify the brightness range by construction, which is the reason they exist. So the workflow is not “threshold the map and see what comes out”, it is:

  1. Threshold the map with something simple and absolute — std below a value you choose, drift below another.

  2. Go class by class, and inside each one draw polygons where the map is quiet.

  3. Check every class still has enough pixels to satisfy whatever minimum your normalization enforces.

Step 1 in code, with thresholds as plain numbers rather than percentiles, so the same criterion means the same thing on another extent, another date or another sensor:

MAX_STD = 0.06      # reflectance units
MAX_DRIFT = 0.02    # reflectance accumulated over the whole series
MIN_OBS = 100       # observations behind the statistics
MIN_PATCH = 25      # pixels; at 10 m that is 2500 m², a 5x5 Sentinel-2 block

quiet = (std.lte(MAX_STD)
         .And(drift.lte(MAX_DRIFT))
         .And(n_obs.gte(MIN_OBS)))

# Isolated pixels are no use as ground targets: they have to survive geolocation
# error and the point spread function of whatever sensor is being normalized
patch_size = quiet.selfMask().connectedPixelCount(maxSize=128, eightConnected=True)
stable = quiet.selfMask().updateMask(patch_size.gte(MIN_PATCH)).rename('stable')

Map = geemap.Map()
Map.centerObject(roi, zoom=8)
Map.addLayer(std, {'min': 0.01, 'max': 0.12,
                   'palette': ['1a9850', 'ffffbf', 'd73027']}, 'std (swir1)')
Map.addLayer(stable, {'palette': ['00d4ff']}, 'quiet and not drifting')
Map

Warning

connectedPixelCount is expensive, and so is anything that forces the whole stack to be evaluated at once. Adding stable to a map is fine — tiles render lazily, and only what you look at gets computed — but a reduceRegion over the full extent on a stack built from eight years of Sentinel-2 will return Computation timed out or User memory limit exceeded, even when every band computes happily on its own. Summarise one layer at a time, or send the job out as an export task.

Both thresholds are yours to set, and the table in section 6 is the calibration: a class you already trust tells you what value your surfaces actually reach. MAX_STD = 0.06 admits everything from open sea to dunes while excluding the gypsum stacks; MAX_DRIFT = 0.02 removes quarries and the stacks while leaving the sea untouched.

7.1 Export the candidates#

Vectorizing an extent this size in the notebook will time out, so it goes out as a task.

polys = stable.reduceToVectors(
    geometry=roi, scale=20, crs=UTM, geometryType='polygon',
    eightConnected=True, labelProperty='stable', maxPixels=1e10,
)

task_vec = ee.batch.Export.table.toDrive(
    collection=polys, description='pia_candidates', folder='ndvi2gif',
    fileFormat='GeoJSON',
)
# task_vec.start()

8. What to do with the output#

The polygons are candidates, not PIAs:

  • Give them a class. The method finds surfaces that do not change; it has no idea whether a given one is a runway, a dune or a factory roof. Normalization schemes that enforce a minimum number of pixels per class need that label.

  • Check the brightness spread survives. Especially after clipping to a sub-area, where a whole class can disappear.

  • Watch for stable-but-wrong surfaces. Deep shadow in relief, permanent cloud over a mountain range and water bodies at a fixed level are all radiometrically quiet and all poor targets, for different reasons.

  • Re-run as the archive grows. Stability is a statement about a time span. A target that was invariant 2018–2025 is not promised to stay that way.

Appendix: selecting by percentile instead#

An alternative to absolute thresholds is to cut the extent into brightness strata and keep the most stable fraction inside each one. It has a real attraction — every stratum contributes candidates, so the brightness spread is guaranteed without any class structure — and a real cost, spelled out at the end.

import numpy as np

BANDS = ['red', 'nir', 'swir1']
N_STRATA, KEEP_PCT, N_SAMPLE = 10, 10, 4000

sigma_b, level_b = {}, {}
for b in BANDS:
    c = series(b)
    sigma_b[b] = c.reduce(ee.Reducer.stdDev())
    level_b[b] = c.mean()

# A PIA has to hold in every band, so the criterion is the worst one
sigma_worst = ee.Image(sigma_b['red']).max(sigma_b['nir']).max(sigma_b['swir1']).rename('sigma')
brightness = level_b['red'].add(level_b['nir']).add(level_b['swir1']).divide(3).rename('brightness')

# A random sample is enough to place the cuts; getInfo() stops at 5000 features
sample = sigma_worst.addBands(brightness).sample(
    region=roi, scale=100, numPixels=N_SAMPLE, seed=42, tileScale=4, dropNulls=True,
).getInfo()
pts = [f['properties'] for f in sample['features']]
sg = np.array([p['sigma'] for p in pts])
br = np.array([p['brightness'] for p in pts])

edges = list(np.percentile(br, np.arange(100 / N_STRATA, 100, 100 / N_STRATA)))
strat = np.digitize(br, edges)
thr = [float(np.percentile(sg[strat == k], KEEP_PCT)) for k in range(N_STRATA)]

stratum = ee.Image(0)
for e in edges:
    stratum = stratum.add(brightness.gt(float(e)))
candidates = sigma_worst.lte(
    stratum.rename('stratum').toInt().remap(list(range(N_STRATA)), thr)
).And(n_obs.gte(MIN_OBS)).rename('candidate')

print('brightness deciles:  ', [round(e, 3) for e in edges])
print('sigma threshold each:', [round(t, 4) for t in thr])

Scored against the same reference classes, with 10 % of the extent selected by construction, and again after also requiring drift below 0.02:

class

% selected

% selected, not drifting

Arena

81 %

62 %

Aeropuertos

76 %

52 %

Urbano-2

66 %

53 %

Urbano-1

58 %

49 %

Marisma mareal

57 %

52 %

Minería

55 %

18 %

Pinar

31 %

23 %

Pastizales

27 %

18 %

Mar

26 %

26 %

Fosfoyesos

19 %

2 %

Embalses

0 %

0 %

Read against the 10 % baseline this is encouraging: runways, urban fabric and dunes come out seven to eight times enriched, and reservoirs are rejected outright — the textbook pseudo-invariant surfaces recovered without being told what they are. The drift column does its job too, collapsing quarries and gypsum stacks while leaving the sea untouched.

And it is the wrong tool for judging a pixel. Every number above is a rank against the rest of this extent. A fifth of it is sea, so the first two brightness deciles are essentially all water and the method ends up sorting calm sea from choppy sea — which is why Mar, the quietest class in the archive, scores a middling 26 %. Clip to a sub-area where the ocean is gone and the cuts move, and the same pixels change verdict without having changed at all.

Percentiles answer “which are the most stable pixels here”. Absolute thresholds answer “has this pixel changed”, which is the question a PIA actually poses, and the only one whose answer travels between extents, dates and sensors. Use the appendix when you want a fixed budget of candidates spread across brightness; use section 7 otherwise.

See also#