Annual cloud-free leaf-on composites of Landsat Collection 2 surface
reflectance for the six New England states and their surroundings, one
composite per year from 2000, as a plain Zarr v3 multiscales pyramid
(30 m to 3840 m, eight levels). Six bands (blue, green, red, nir, swir1, swir2) as uint16 digital
numbers with reflectance = DN * 0.0000275 - 0.2, plus clear_count, the
number of clear looks each pixel's composite drew on, and three provenance
planes (source, pick_doy, borrowed_pct) that say, in the tile-years
where the build lent looks from the neighbouring years, which year and
day each pixel came from.
Grid: EPSG:4326, 0.00025 degree pixels, 26200 x 27600, bbox -73.75, 40.95, -66.85, 47.5. Window: June 1 to September 30 where a pixel has at least four clear looks, otherwise May 1 to October 31.
Data: CC0-1.0. Code: MIT. Source: Landsat Collection 2 Level-2 Tier 1 via the Microsoft Planetary Computer. The code that built it, with the measurements behind every rule here, is at https://github.com/kentstephen/ne-landsat-temporal-composite.
The Source Cooperative repository ID is
kentstephen/landsat-mosaics-new-england. This project set out to be
a mosaic. It is now a composite, not a mosaic.
A mosaic joins scenes that cover different ground; a composite reduces
scenes that cover the same ground on different dates to one value per
pixel, which is what the medoid does here. Read "mosaic" in the ID, and
in the older code and notebook names, as "composite".
Everything is public and read over https, no credentials. Set the prefix once:
BASE = "https://data.source.coop/kentstephen/landsat-mosaics-new-england"
The pyramid with zarr-python. Consolidated metadata sits at the root, so open the root and index into a level. Level 0 is 30 m; level 7 is 3840 m and fits in memory for the whole region.
The pyramid with xarray. One level is one dataset with time, y,
x coordinates. Leave consolidated at its default (or True): plain
https cannot list a directory, so the level is found through the
consolidated metadata at the root, and consolidated=False opens an
empty dataset.
Zeros in the bands are nodata (fill 0), and borrowed_pct == 255 means
outside the ladder tile-years. Mask both before averaging.
The Icechunk store. build.icechunk/ is level 0 only, opened by tag.
The proxy speaks S3 with the account as the bucket, anonymous, path
style.
To see what the ladder changed, read the same window from v3 (before)
and v4 (after) and diff where source is 1 or 2.
Currently 2000 to 2025, 26 years. A new year is appended to this same
store once its leaf-on window is fully acquired, so the length of the
time axis grows; read it from the store rather than assuming a count.
Appends add the new year only and never change published pixels. If the
method changes (reducer, mask, window, fill rule, grid, sensors), the
result is published as a new pyramid with a new name rather than an
append here. Each append is listed below with its date and the Icechunk
tag it was exported from, so a reader can say which state they used.
2000 to 2025 initial release tag v4-provenance
Read these before comparing years or measuring anything.
One reducer, medoid, every year. Each pixel is one real
observation: the clear look nearest the per-pixel median in band space,
so a pixel's spectrum is never a blend of dates. The store attribute
reducer_by_year records this per year.
Landsat 7 scan line corrector gaps, 2003 to 2023. From June 2003
every Landsat 7 scene has wedge-shaped gaps about 1 km apart. Inside a
gap zone a pixel has a different set of clear looks than its
neighbours, and any per-pixel reducer picks different dates in adjacent
stripes, which showed as diagonal tone striping in every year with
Landsat 7 looks (2003 to 2023), strongest 2003 to 2012. The build
therefore limits Landsat 7 to a fill role in those years, with one of
two rules per tile-year (4096 px tiles, 693 tile-years in 2003 to
2023). The fill rule, 314 tile-years: at each pixel the Landsat 7 looks
count only where Landsat 5, 8 and 9 together give fewer than four clear
looks. The ladder, the 379 tile-years that still striped under the
fill rule: own-year Landsat 5, 8 and 9 looks where they reach three;
else those plus the same platforms' looks from the year before and
after, where together they reach three; else every own-year look,
Landsat 7 included, plus the neighbour-year looks. Where the other
platforms suffice, the gap pattern does not change the look set and
the striping goes; where they do not, Landsat 7 fills in and the count
still shows the gap geometry. The store attributes l7_fill_by_year,
min_clear, l7_ladder_tile_years and l7_ladder_note record which
rule applies where. 2000 to 2002 (Landsat 7 before the failure) and
2024 onward (no Landsat 7) were built without either and have no
SLC-off looks. Measured with stripe_score.py (directional power at
the gap period, NIR band, per 30 km block of level 1, median over
blocks; 1.0 is isotropic, years with no Landsat 7 sit at 1.05 to 1.10):
The cost is fewer looks per pixel where Landsat 7 was excluded (mean clear count 8 to 10 down to 5 to 6 in 2003 to 2021), and in ladder tile-years a pixel may rest on three looks rather than four.
In 379 tile-years a pixel may be an observation from the year before or after. The ladder's second rung lends neighbour-year looks, so in those tile-years the composite year is not guaranteed per pixel: a change on the ground can show a year early or late there. Each pixel is still one real observation, never a blend. Over the pixels in the 379 tile-years, 92.9% are own-year Landsat 5, 8 or 9 picks, 2.6% the year before, 3.4% the year after, 0.7% own-year Landsat 7. Per year the borrowed share runs from 0.0% (2012, 2023) to 14.9% (2013, all from 2014, since 2012 has nothing to lend) and 14.7% (2019). The planes say where: (level 0) classes each pixel as own-year, year before, year after or own-year Landsat 7; gives the day of year of the stored observation, in that observation's own year; (every level) is the share of a pixel's footprint that came from a neighbour year, 0 or 100 at level 0, and 255 outside the ladder tile-years. A hold-out test on seven windows () measured each rung against a well-sampled own-year medoid, median |NDVI error|, split by whether the landscape changed between the neighbour years:
Chunks are 1 x 512 x 512 at every level, in 1 x 4096 x 4096 shards on levels 0 to 2, zstd. Consolidated metadata is at the root.
Beside the pyramid sits build.icechunk/, the Icechunk store the pyramid
was exported from: level 0 only, the same arrays (no borrowed_pct),
readable with icechunk-python. It is the working store for the yearly
append, not the product, and it holds three tags:
v3 is kept so the pre-ladder pixels of the 379 tile-years stay readable;
source says which pixels were borrowed but not what was there before.
Nothing older is kept.
Three vector layers ship alongside the store in supplemental/, each
with its source, license and transformations in supplemental/README.md:
The composites are built from USGS Landsat Collection 2 Level-2 Tier 1 surface reflectance, read through the Microsoft Planetary Computer.
tiger_states.parquet is from the U.S. Census Bureau's TIGER/Line 2024
state shapefile, and nhd_water_bodies.parquet from the USGS National Hydrography
Dataset High Resolution (December 2023), both public domain.
hf435_wildlands.parquet is redistributed from the Harvard Forest Data Archive,
reprojected and simplified. Credit the original:
Wildlands of New England GIS Data 1900-2022. License CC0 1.0. Foster D, Johnson E, Hall B. 2023. Wildlands of New England GIS Data 1900-2022. Harvard Forest Data Archive: HF435. https://harvardforest1.fas.harvard.edu/exist/apps/datasets/showData.html?id=HF435 Compiled by Wildlands, Woodlands, Farmlands & Communities with Harvard Forest, the Highstead Foundation and the Northeast Wilderness Trust.
sourcepick_doyborrowed_pctstats/holdout/Landsat 5 scene edges are eroded. From 2002 Landsat 5's scan mirror ran in bumper mode and its east and west scene edges are combs, one tooth per 16-line swath, kilometres long. In path overlaps the look set alternated swath by swath and the medoid striped with it. Each Landsat 5 look's data-present mask is eroded by 25 grid rows (about 700 m) before compositing, so pixels near those edges draw on the neighbouring path.
Landsat 7 and 9 same-day looks, 2022 to 2023. Landsat 7's orbit was lowered in 2022 and it drifted to a pass that images a tile on the same day as Landsat 9 (60 tile-years). The two are separate looks, keyed by day and platform, so the Landsat 7 gaps stay under the fill rule.
2012 draws on 2011 and 2013. Landsat 5 had stopped and Landsat 8
had not launched, so 2012 has Landsat 7 alone and the fill rule would
have nothing to fall back on. Its composite pools the 2011 to 2013
clear looks under the rule, so a 2012 pixel may be a 2011 or 2013
observation. The store attribute year_span_by_year records this.
Pooled, 2012 has a mean of about ten clear looks per pixel and scores
like a year with no Landsat 7; on its own looks it had five and was the
most striped year.
2003 has holes. Where every Landsat 7 look fell in a gap and every
Landsat 5 look was cloudy, the 2003 composite has no data (all bands 0,
clear_count 0): about 0.05% of the area, in runs up to about 1 km
wide. The fill rule does not change this (no other look existed). No
other year has inland nodata.
30 m source on a finer grid, resampled nearest. The 0.00025 degree pixel is 27.8 m north-south and about 20 m east-west at 43.5N, so each 30 m source pixel covers about 1.1 by 1.5 grid pixels and roughly one column in three duplicates its neighbour. There is no resolution beyond 30 m in the data. For a true 30 m projected grid, reproject with an averaging resampler.
Water reflectance goes negative. Collection 2 aerosol correction over water and dark shadow produces reflectance below 0.0 (DN below 7273), mostly in nir and swir: about 3% of valid nir pixels, nearly all over the ocean and lakes. This is the USGS product, not a compositing artifact, and it is kept as is.
Saturation. A few hundred thousand pixels across all years hold 65535 in one band (bright targets, most in 2012 and 2003). Values above 65455 are outside the Collection 2 valid range.
Water is not masked. The composites cover the full bbox including the Gulf of Maine, Long Island Sound, and adjacent Quebec, New Brunswick, and New York.
Pyramid levels are means. Level N is the nodata-aware mean of
2^N by 2^N level 0 pixels, rounded half up, in the same DN scale. Coarse
levels of 2003 look filled where level 0 has holes, because the mean
bridges them. clear_count is averaged the same way.