Lab 4: Building and Reducing Image Collections
Filter, map, and reduce: from raw archive to composite image
Background
Every Landsat image in Google Earth Engine is a single moment — one overpass, one cloud geometry, one atmospheric state. In Labs 2 and 3 you worked with single scenes, and that was appropriate: both labs were about reading what a sensor sees at a specific time and place. But most environmental questions are not about a moment. They are about a process. A drought unfolds over months. A forest is cleared across multiple growing seasons. A river delta shifts through years of sediment delivery. Answering those questions requires not one image but a collection of images — and a principled way to compress that collection into something analytically useful.
The Landsat archive holds more than 1.9 million scenes from the Landsat 8 instrument alone. That number is not useful in its raw form; it needs to be filtered, preprocessed, and then statistically summarized. Cardille et al. (Ch. F12) frames this as a three-step workflow: filter the archive to a manageable subset, map a preprocessing function over every image in that subset, and reduce the result to a composite that represents what the scene typically looks like. Each step is not arbitrary software design — each solves a specific scientific problem that arises when you try to turn a pile of imagery into an interpretable result.
Before any pixel value can be used in a statistic, it must mean something in physical units. Landsat Collection 2 Level 2 products store surface reflectance as 16-bit unsigned integers — a compact encoding that must be converted to true reflectance values before any comparison across dates or locations is valid. Campbell Ch. 11.4 describes this conversion as a mandatory radiometric preprocessing step; the scale factor is not optional decoration but a prerequisite for any meaningful calculation. And before averaging across a stack of images, pixels contaminated by clouds must be identified and excluded. A single cloud in fifteen scenes can double the apparent reflectance of the pixel it covers, pulling statistics away from the surface they are supposed to describe. The pixel-by-pixel quality flags bundled with each Landsat acquisition identify cloudy pixels before they can corrupt the result.
By the end of this lab you will have built a complete processing chain from the raw archive to a monthly NDVI time series, applied it to two locations with very different phenological stories — Lisbon, Portugal, and the Skagit Valley in Washington — and then applied the same chain independently to a third location, the Gran Chaco of Paraguay. That last application asks you to interpret a record from 2020, a year when fire and drought produced a signal that stands apart from the seasonal vegetation cycle and connects directly to the deforestation history Hansen et al. (2013) documented from space.
Part 1 — Building and Exploring a Collection
Filtering the archive
The Landsat 8 Collection 2 Level 2 archive contains nearly two million images acquired since the sensor’s launch in 2013. That full archive is the starting point for every analysis — but it is also useless until it is narrowed down to the images that are relevant to a specific question. Script L4.1 makes the narrowing process explicit by printing the size of the collection after each filter is applied. You will watch the count fall in three stages: from millions down to a few hundred (by location), then to a few dozen (by year), then to a working set of about fifteen (by scene-level cloud cover).
The order in which those filters are applied is not arbitrary. Cardille et al. (Ch. F12) explain the efficiency argument: filtering by location first eliminates the overwhelming majority of images before the date filter even runs, because Landsat 8 has imaged most of the Earth’s surface thousands of times but each individual tile has been imaged roughly once every sixteen days. Applying the most selective filter first is the difference between a query that returns in seconds and one that times out. The script demonstrates both orderings — efficient and inefficient — and prints identical final counts for both. That identical count is the confirmation that filter order changes processing speed, not the result.
The scene-level cloud cover threshold you apply in this step does something conceptually different from the pixel-level masking you will do in Part 2. It operates on a single number attached to each scene as metadata — an estimate, made by the USGS ground processing system, of what fraction of the scene is covered by cloud. Filtering on this number removes entire acquisitions from the collection; it says nothing about which individual pixels are cloudy. A scene that passes a 50% threshold may still have the exact location you care about sitting under a cloud. That is why the scene filter and the pixel mask are both necessary and why neither replaces the other.
Visualizing what is in the collection
Script L4.2 builds on the filtered collection from L4.1 and introduces two ways to look at it before reducing. The first is a per-pixel count map — a spatial image showing how many times each geographic location was captured across all the scenes that passed your filters. The count is not uniform; it depends on orbit geometry, the overlap between adjacent Landsat swaths, and how many of the scenes over a given pixel passed the cloud threshold. Locations near the edge of the swath, where adjacent orbits overlap, tend to have higher counts. Densely cloudy regions over the ocean or near coastlines tend to have lower counts.
The second is a time series chart — a plot of one band’s value at a specific point across all the acquisition dates in the collection. This chart is not yet a clean phenological record; it samples raw, unscaled digital numbers that include both the surface and any cloud contamination that slipped past the scene-level filter. The variation in the chart is diagnostic: dates where the value spikes sharply above the surrounding observations correspond to acquisitions where the sample pixel was covered or partially covered by cloud. Campbell Ch. 11.2’s description of the histogram as a tool for diagnosing value distributions applies directly here — the time series gives you the same distributional information but spread across a temporal rather than spatial axis.
Part 2 — Preprocessing Every Image
Why the values must change before you reduce
The pixel values stored in a raw Landsat Collection 2 Level 2 image are integers in the range of approximately 7,000 to 60,000. Those numbers do not represent reflectance; they represent reflectance encoded as integers through a standard affine transform. The true surface reflectance — a dimensionless ratio between zero and one — can only be recovered by applying the official USGS scale factor: multiply by 0.0000275 and subtract 0.2. This is the conversion Campbell Ch. 11.4 describes as radiometric preprocessing. Until it is applied, the pixel values are internally consistent but physically uninterpretable: an average of raw integers is a mathematical result, not a meaningful estimate of surface conditions.
Cloud masking operates on a different band entirely. The Landsat C2 L2 product includes a QA_PIXEL band alongside the spectral bands. Each pixel in QA_PIXEL is a 16-bit integer whose individual bits encode separate quality conditions: whether the pixel is cloud, cloud shadow, dilated cloud, or cirrus. The mask is constructed by reading those bits and zeroing out pixels where any of the relevant flags are set. This is a per-pixel operation — not a scene-level judgment, but an image-by-image, pixel-by-pixel quality screen applied at the full spatial resolution of the data. Pixels that fail the mask become transparent and are excluded from any subsequent statistics.
Mapping a function over the collection
Script L4.3 introduces the map step of the filter-map-reduce workflow. Rather than writing a loop that processes one image at a time, a single preprocessing function — maskAndScale() — is defined and then applied to every image in the collection in one line. The result is a new collection of the same size, where every image has been cloud-masked and radiometrically scaled. The .map() pattern is not just a programming convenience; it is the foundation of Earth Engine’s distributed computing model. Earth Engine parallelizes the map step automatically across its compute infrastructure, so a collection of fifteen images is processed no differently from a collection of ten thousand.
The before-and-after map display in this script lets you toggle between the raw, unscaled collection mosaic and the processed collection. The processed layer will have transparent gaps where clouds have been removed — the Google satellite basemap shows through those gaps. The values visible in the processed layer are in the range of 0.0 to about 0.4 for typical land surfaces, which is the reflectance range Campbell Ch. 11.4 describes for vegetated and urban cover.
Part 3 — Reducing a Collection to a Composite
Three reducers, three assumptions
Once the collection has been preprocessed, it can be reduced — collapsed from a stack of fifteen images into a single composite image where each pixel holds one representative value. Script L4.4 produces three composites using three different reducers, and the differences between them reveal something important about the distribution of values in the stack.
The mean composite assigns equal weight to every observation. In a cloud-free stack, this is a sensible and statistically efficient choice. In a stack where some cloud contamination has survived the mask — thin cirrus, cloud edges, shadow boundaries where the QA flag did not fire — the mean is the most vulnerable. A single anomalously bright cloud-edge pixel in fifteen observations pulls the mean upward without any visible warning in the output.
The median composite picks the middle value when all observations at a pixel are ranked in order. It is robust to outliers by construction: whether the contaminating value is the highest in the stack or the second-highest, it does not affect the median as long as fewer than half the observations are contaminated. For the fifteen-image Lisbon collection, median and mean often agree closely — but where they diverge, the gap is diagnostic of contamination that survived the mask.
The 30th-percentile composite — or whichever percentile you choose in the script’s CONFIGURATION — deliberately picks a value from the lower end of the ranked distribution. The logic here, which Cardille et al. (Ch. F13) explain directly, is asymmetric: cloud contamination always pushes pixel values upward (clouds are bright), never downward. A percentile from the lower portion of the distribution is therefore systematically less likely to include contamination. The tradeoff is that it also picks genuinely lower-reflectance observations — dark atmospheric conditions, shadows, or partially senesced vegetation — over brighter, greener ones.
Reading the histogram
The script’s final section places all three composite distributions side by side in a histogram of the red band (SR_B4) over a 50-kilometer region around the Lisbon study point. The histogram makes the distributional logic visible rather than abstract. Before reading any explanations of what you are supposed to see, examine the histogram and describe what you observe: where does each distribution sit relative to the others, and what is the shape of each curve?
Part 4 — Monthly Aggregation and Vegetation Phenology
From one composite to twelve
A single annual composite collapses an entire year of observations into one image. That is useful for mapping land cover, but it loses the temporal signal — the seasonal rhythm of greening and senescence that distinguishes a wheat field from a conifer forest, or a wet-season tropical forest from a drought-stressed one. Temporal aggregation recovers that signal by dividing the year into windows and reducing each window independently.
Script L4.5 moves from Lisbon to the Skagit Valley in Washington, a broad agricultural basin flanked by the Cascade Range. The script partitions the 2022 calendar year into twelve one-month windows and applies the same preprocessing and median reduction to each window. The result is a collection of twelve monthly composites, each tagged with its acquisition month so that a time-series chart can place them correctly on the time axis.
The NDVI — Normalized Difference Vegetation Index — is computed for each image before the monthly reduction. NDVI measures the contrast between red and near-infrared reflectance; live green vegetation absorbs strongly in the red and reflects strongly in the near-infrared, producing values close to 1. Bare soil, water, and impervious surfaces produce lower values; dense green canopy approaches 0.9. Campbell Ch. 17.3 describes vegetation phenology as one of the primary applications of multi-temporal NDVI data: the annual curve of NDVI at a single point encodes the timing of green-up, peak biomass, and senescence in a form directly readable from the chart.
Interpreting the Skagit Valley chart
The NDVI time series you produce in this script reflects conditions at a single point — a specific pixel in the Skagit Valley lowlands. The shape of that curve is determined by what the pixel overlaps, and in an agricultural valley with cover crops, row crops, fallow fields, and managed wetlands, different pixels produce very different annual rhythms. Before drawing conclusions about regional phenology from a point chart, it is worth examining the NDVI spatial composite for the same months and checking whether the sampled pixel is representative of the surrounding landscape or an outlier.
Not every month will produce a visible data point on the chart. The Pacific Northwest receives substantial cloud cover through fall and winter; some monthly windows may contain no clear-sky observations at the sample point after masking, and those months will appear as gaps. This is not an error — it is an accurate representation of data availability, and it is the same information limitation that Hansen et al. (2013) had to manage when constructing their global forest loss product from Landsat data.
Part 5 — Independent Application
The Gran Chaco, 2020
Script L4-EX1 asks you to apply the complete filter-map-reduce-aggregate pipeline — the same engine you used in Parts 1 through 4 — to a pre-assigned location and year. The study area is in the Gran Chaco of Paraguay, one of the four global deforestation hotspots documented by Hansen et al. (2013) in the paper you read for this lab. The year is 2020.
The Gran Chaco is a dry subtropical forest with a pronounced wet season from November through April and a dry season from May through October. At this location, the landscape is a mosaic of intact forest and converted agricultural land — the geometric pattern of clearing is visible in the satellite basemap and in the NDVI composite you will produce. Hansen et al. used Landsat imagery and the same filter-map-reduce logic underlying these scripts to quantify how rapidly that mosaic was shifting toward agricultural land between 2000 and 2012.
The year 2020 adds a second layer to the analysis. In that year, the Chaco and the adjacent Pantanal wetland experienced one of the most severe droughts and fire seasons on record. The NDVI time series you produce should show two things: the normal seasonal arc driven by the wet-dry cycle, and an anomalous departure in one or more of the late dry-season months. Your task is to describe what you see, identify where the anomaly falls in the seasonal cycle, and connect it to what you know about NDVI and vegetation stress.
Your two remaining design decisions — the scene-level cloud threshold and the peak-growing-season month — are yours to choose and justify. The script’s guidance on both reflects what you have learned about the relationship between scene-level filtering and pixel-level masking, and about the southern-hemisphere timing of peak canopy cover. Read that guidance before choosing your values, and justify your choices in your SOP.
Interpretation Questions
Your answers to these questions belong in the Results section of your SOP.
Q1 (Script L4.1) Campbell Ch. 11.2 defines the mean of a set of values as the sum divided by n, where n is the number of observations. In the filter-map-reduce workflow, the filter chain determines what n will be for every statistic computed later. Using your Console output from Script L4.1, describe how n changed at each of the three filter stages. Then, using the milk-store analogy Cardille et al. (Ch. F12) introduce, explain what each of your three filters corresponds to in that analogy. What are the “stores,” and what does each filter represent as a step in narrowing the list?
Q2 (Script L4.2) Campbell Ch. 11.2 introduces the histogram as a chart of how pixel values are distributed within a single band. The time series chart you produced in Script L4.2 is also a distribution — but organized along a different axis. Describe what axis a histogram uses and what axis the time series chart uses. Then examine your time series chart and identify at least one acquisition date where you can tell from the chart alone that cloud contamination likely affected the sample pixel. Explain what feature of that data point tells you the pixel was contaminated, using the reflectance properties of clouds that Campbell Ch. 2 describes.
Q3 (Script L4.3) Campbell Ch. 11.4 describes the conversion from digital number to surface reflectance as a prerequisite for any physical interpretation of pixel values. Explain in your own words why applying
ee.Reducer.mean()to the raw, unscaled collection — the BEFORE layer in Script L4.3 — would produce a result that cannot be compared to the AFTER composite, even if both images look similar on the map. In the same answer, explain why cloud masking is not optional before reducing: use a specific value from your time series chart in Q2 to illustrate what would happen to the mean at that pixel if the contaminated observation were included.
Q4 (Script L4.4) Examine the three SR_B4 values printed to the Console in Script L4.4 — the mean, median, and 30th-percentile composites at the Lisbon point. Campbell Ch. 11.2 distinguishes the mean from the median by their sensitivity to outliers. Using those three values, describe (a) what the gap between mean and median tells you about the distribution of pixel values in the stack at this location, and (b) why the 30th-percentile value is systematically lower than both. Then look at the histogram the script produces. Use the legend colors to identify each curve and describe the shape and relative position of the three distributions: where does each distribution peak along the reflectance axis, which curve sits farthest to the left (toward lower reflectance), and how do their spreads compare? Connect what you observe to the compositing logic that produced each reducer — why would one distribution be systematically shifted toward lower values?
Q5 (Script L4.4) Script L4.4 includes a parameter —
CLOUD_MAX_SCENE— that screens entire Landsat scenes before any pixel-level processing occurs, and a bitmask insidemaskAndScale()that removes individual cloud pixels. These are two distinct quality controls operating at two different scales. Explain the difference in what each control does and what each one cannot do alone. At what point does changingCLOUD_MAX_SCENEaffect the chart at a single pixel, and when does it not? Your answer should draw on what you observed when experimenting with the threshold in Script L4.4.
Q6 (Script L4.5) The NDVI time series chart from Script L4.5 samples a single pixel in the Skagit Valley. Describe the shape of the annual NDVI curve you see: when does the vegetation at your sample point reach its highest values, when does it reach its lowest, and how abrupt or gradual are those transitions? Campbell Ch. 17.3 discusses phenological cycles in terms of the underlying biological and climatic drivers. Based on those drivers, suggest what land cover type or crop cycle is most consistent with the timing and amplitude of the curve you observe. Then compare the timing of peak NDVI in your Skagit Valley chart to the timing of peak NDVI in the Gran Chaco chart from L4-EX1 — what does the difference tell you about the two regions?
Q7 (Script L4-EX1) Examine the peak-month NDVI map and the monthly time series chart from Script L4-EX1. In the map, describe the spatial pattern of NDVI values — where are the highest values, where are the lowest, and what land cover types correspond to each when you toggle to the satellite basemap? Draw on the spectral vocabulary you developed in Module 3 when describing what you see. In the time series chart, identify any month where NDVI departs from the seasonal arc, and describe the departure in terms of timing, magnitude, and what it suggests about the nature of the disturbance — vegetation stress versus vegetation removal. Explain the most likely physical cause using what you know about how drought and fire affect surface reflectance, and connect both observations to what Hansen et al. (2013) documented about forest loss in this region.
Standard Operating Procedure
Complete this section in your own words as you work through the lab. Your SOP documents your workflow so that a technically competent colleague could replicate your exact analysis.
Study Area
This lab spans three geographic locations used across the five scripts. Identify each study area by name, coordinates, and a brief description of the environmental or land use context that makes it relevant to the questions this lab asks.
Data Sources
List every dataset loaded across the five scripts. For each entry, record the dataset name, collection ID as it appears in GEE, sensor, spatial resolution, and the date range you filtered to.
| Dataset | Collection ID | Sensor | Spatial Res. (m) | Date Range |
|---|---|---|---|---|
Add rows as needed.
Methods
Describe the workflow you applied in each script. Your Methods section should explain how you built the filter chain, what function you mapped over the collection and why, which reducer you chose for each composite and the reasoning behind it, and how you structured the monthly aggregation. For the independent exercise, explain the design decisions you made for cloud threshold and peak month, connecting each choice to the regional climate and sensor logic described in the script.
Results
Include the following in this section:
- Screenshots of the BEFORE and AFTER mosaics from Script L4.3, with descriptive captions
- The Console output showing mean, median, and 30th-percentile SR_B4 values at the Lisbon point from Script L4.4
- The SR_B4 histogram comparing all three reducer distributions from Script L4.4
- The NDVI time series chart from Script L4.5 (Skagit Valley)
- The peak-month NDVI map and NDVI time series chart from Script L4-EX1 (Gran Chaco)
- Answers to Interpretation Questions Q1–Q7
Limitations
Identify the limitations of the analysis you performed in this lab. These should be methodological or data-related — not personal difficulties with the software.
Consider: the spatial scale at which a single-pixel chart represents a landscape; what the scene-level cloud threshold can and cannot protect against; the choice of calendar year for monthly aggregation and what it excludes; whether a single-point NDVI sample is representative of the broader land cover visible in the map; and what would be needed to distinguish fire, drought, and land use change as causes of an NDVI anomaly.
Submission Checklist
ESCI 442/542 — Introduction to Remote Sensing | Western Washington University