Assignment 6#

Due: Wed Oct 28th at 11:59 pm ET

The goal of this assignment is to retrieve and analyze a seven-year time series from two very different kinds of big Earth data, and to do it within a fixed memory budget. You will build a satellite time series of reservoir water area from Sentinel-2 imagery with stackstac and Dask Array, retrieve daily precipitation records from NOAA’s GHCN-D archive with Dask DataFrame, and put the two together to ask a simple question: how did the surface area of Folsom Lake, California, change between 2019 and 2025, and how does that track with local precipitation?

You will write the imagery pipeline by hand, then use an AI coding agent to build the precipitation half. As in the previous assignments, your main job on the agent’s work is to catch where it is wrong, and in this assignment “wrong” includes code that runs out of memory or quietly pulls far more data than it needs.

What you are practicing

  • Dask Array through stackstac: lazy stacks, chunks, and reducing before you compute.

  • STAC time-series retrieval: filtering a search, removing duplicate items, and requesting only the assets, resolution, and data type you need.

  • Dask DataFrame: lazy Parquet reads from S3, column selection, early filtering, persist, and repartition.

  • Memory discipline: estimating the size of a computation before running it, and measuring what it actually used.

  • Judgment: reviewing agent-written code for silent data errors and memory waste.

Notebooks are for visualization, not logic

As in the previous assignments, your notebook must not contain your analysis functions. Every function you write, and every function the agent writes, lives in src/utils.py. The notebook imports those functions and uses them to explore, measure, and plot.

All work lives in a new directory named assignment-6/ in your geog313-assignments repository. Follow the same git discipline as before: stage files explicitly by name (never git add .), commit in small focused steps with clear messages, and push as you go. Read the whole assignment before you start, including the background on memory below.

Background: why Dask workflows run out of memory#

Dask is lazy, so building a stack of 1,000 satellite scenes costs almost nothing. Memory problems show up only when the computation runs, and in practice they almost always come from one of three places.

The first is data that is bigger than it needs to be before any computation starts. The size of a stackstac array is simply

n_time x n_band x n_y x n_x x bytes_per_value

and every term in that product is under your control. Duplicate items inflate n_time. Stacking every asset instead of the two or three you need inflates n_band. Reading at 10 m when 20 m is enough multiplies n_y x n_x by four. The stackstac default data type is float64 (8 bytes per value), while float32 (4 bytes) is more than enough precision for reflectance. These choices multiply, so together they can easily make the same analysis a hundred times larger.

The second is bringing unreduced data back to your notebook. Calling .compute(), .values, .load(), or .plot() on a full stack asks Dask to assemble every pixel of every scene in the memory of your Python kernel. The pattern that scales is the opposite: reduce lazily (a mean, a sum, a count over x and y) and compute only the small result.

The third is the cluster itself. Each time you run a cell that creates a LocalCluster, you start a new cluster with its own workers, and the old one keeps its memory until the kernel dies. A few re-runs on an 8 GB laptop is enough to exhaust it.

Duplicate items in the Planetary Computer catalog#

Two features of the Planetary Computer sentinel-2-l2a collection make the first problem worse over our area of interest:

  • Sentinel-2 is distributed in 110 km MGRS tiles that overlap their neighbors. Folsom Lake sits near the northern edge of tile 10SFH, inside its overlap with tile 10SFJ, so every acquisition is returned twice, once per tile.

  • Some scenes have been reprocessed by ESA with a newer processing baseline, and the catalog keeps both versions. For most dates between February 2022 and December 2023 you will find, for each tile, the original item (baseline 04.00 or 05.09) and a reprocessed one (baseline 05.10) with the same acquisition time.

As a result, a search over this AOI returns about two-thirds more items for 2023 than for 2024, even though the satellites observed the lake on the same schedule. Every item property you need to sort this out (s2:mgrs_tile, s2:processing_baseline, eo:cloud_cover) is in the item metadata; look at it before you stack anything.

From digital numbers to reflectance#

Sentinel-2 L2A pixels are stored as integer digital numbers (DN). Since processing baseline 04.00 (introduced on January 25th, 2022), ESA adds an offset of 1000 to every reflectance DN, so that

reflectance = (DN - 1000) / 10000    for baseline >= 04.00
reflectance = DN / 10000             for baseline <  04.00

The Planetary Computer items do not carry this scale and offset in their metadata, so stackstac cannot apply it for you; you have to apply it yourself, scene by scene, based on each item’s s2:processing_baseline. Do not skip this because NDWI is a ratio: the ratio cancels the scale, but not the offset. Ignoring the offset shifts NDWI for every scene from 2022 on, and subtracting it from a scene that does not have it pushes dark water below zero reflectance and flips the sign of its NDWI. Your series will still look plausible at a glance.

Data and area of interest#

Area of interest (AOI), given as a bounding box (min_lon, min_lat, max_lon, max_lat) in EPSG:4326:

(-121.18, 38.69, -121.04, 38.81)

This box covers Folsom Lake and lies entirely inside MGRS tile 10SFH, which is in UTM zone 10N (EPSG:32610). Use the period January 1st, 2019 to December 31st, 2025 for the imagery. Folsom Lake’s surface area at full pool is about 11,450 acres (about 46 km²), which you will use as a sanity check.

Imagery comes from the Microsoft Planetary Computer STAC API (https://planetarycomputer.microsoft.com/api/stac/v1), collection sentinel-2-l2a, as in Lecture 13. Its assets are named by band number: green is B03, near infrared is B08, and the Scene Classification Layer is SCL.

Precipitation comes from the NOAA Global Historical Climatology Network Daily (GHCN-D) archive on AWS (s3://noaa-ghcn-pds, public, no account needed), which you used in Lecture 14. This time, use the Parquet copy rather than the CSVs, which is partitioned by year and by element:

s3://noaa-ghcn-pds/parquet/by_year/YEAR=2023/ELEMENT=PRCP/*.parquet

Station locations are in https://noaa-ghcn-pds.s3.amazonaws.com/ghcnd-stations.txt (a fixed-width text file), and the format, units, and quality flags of every variable are documented in https://noaa-ghcn-pds.s3.amazonaws.com/readme.txt. Read the readme before you write any code that touches the values.

Project layout#

Your assignment-6/ directory must follow this layout:

assignment-6/
  pixi.toml            environment, dependencies, and tasks
  pixi.lock            locked versions (managed by pixi)
  .gitignore
  src/
    __init__.py        makes src an importable package
    utils.py           all analysis functions
  notebooks/
    analysis.ipynb     imports from src.utils; explores, measures, and plots only
  tests/
    test_utils.py      tests for the functions in src/utils.py
  results/             exported tables and figures
  AGENTS.md            guidance for the agent
  README.md            what this is and how to reproduce it
  reflection.md        your written reflection (no AI)

The memory budget#

Your pipeline must run within 4 GiB

Start exactly one Dask cluster per notebook, at the top, with this configuration:

from dask.distributed import Client, LocalCluster

cluster = LocalCluster(n_workers=2, threads_per_worker=2, memory_limit="2GiB")
client = Client(cluster)
client

This is deliberately smaller than most laptops, so that the assignment runs the same way for everyone. Run this cell once. If you need a fresh cluster, restart the kernel rather than re-running the cell, and end your notebook with client.close() and cluster.close(). Keep the Dask dashboard open while you work; the memory bars per worker will tell you immediately when something is wrong.

If a worker is paused, restarted, or killed (you will see warnings such as “Worker is at 80% memory usage” or a KilledWorker error), do not raise the memory limit. Go back to the size formula above and find which term is larger than it needs to be, or which cell is computing something unreduced.

Part 0: Reproducible environment (5 pts)#

  1. From the root of your geog313-assignments repository, create the project assignment-6 and initialize it with pixi.

  2. Set the platforms in pixi.toml to ["osx-arm64", "osx-64", "linux-64"].

  3. Add these packages: pystac-client, planetary-computer, stackstac, xarray, dask, distributed, bokeh, pandas, pyarrow, s3fs, pyproj, matplotlib, jupyterlab.

  4. Create the src/, notebooks/, tests/, and results/ structure above, add src/__init__.py, and write a README.md that says what the project does and gives the command(s) to reproduce it.

  5. Commit the scaffolding.

Part 1: The imagery time series, by hand (40 pts)#

Do not use the agent for this part. You are strictly prohibited from using agentic coding here. Writing the imagery pipeline yourself is what lets you judge the agent’s Dask code in Part 2.

In src/utils.py, write and document these functions. All of them must stay lazy: none of them may call .compute(), .values, .load(), or .persist().

  1. search_items(bbox, start, end, max_cloud): search the Planetary Computer with pystac-client for sentinel-2-l2a items that intersect bbox between start and end, keeping only items with eo:cloud_cover below max_cloud. Do the cloud filtering in the search itself (with the query argument of Client.search), not after the items have been downloaded. Return the item collection.

  2. dedupe_items(items, mgrs_tile): keep only items from the given MGRS tile (for example "10SFH", matched against s2:mgrs_tile), and for each acquisition date keep only the item with the highest s2:processing_baseline. The baseline is stored as a string such as "05.10"; compare it as a number. Return the reduced list of items.

  3. to_reflectance(stack): convert a stack of DNs to surface reflectance, scene by scene, as described in the background section, leaving the scl band unchanged. stackstac exposes each item’s s2:processing_baseline as a coordinate along time; use it rather than a hard-coded date.

  4. stack_items(items, assets, bbox, resolution=20, epsg=32610): build a lazy reflectance stack with stackstac.stack that contains only the requested assets, is clipped to bbox with bounds_latlon, uses the given resolution and epsg, and uses dtype="float32" with a float32 NaN fill_value. Pass rescale=False: the items have no scale or offset metadata to apply, and stackstac refuses rescale=True together with float32. Take assets as common band names ("green", "nir", "scl"), map them to the asset keys (B03, B08, SCL), label the band dimension with the common names, and return the result of to_reflectance.

  5. mask_clouds(stack, bad_classes=(0, 1, 3, 8, 9, 10)): use the Scene Classification Layer (scl) to set every pixel whose class is in bad_classes to NaN in the remaining bands. These classes are no data, saturated or defective, cloud shadow, cloud medium probability, cloud high probability, and thin cirrus (see the SCL table). Return the masked stack without the scl band.

  6. ndwi(stack): compute the normalized difference water index, (green - nir) / (green + nir), as a labeled DataArray.

  7. water_area(ndwi, threshold=0.0): for each time step, return a lazy Dataset with two variables: water_km2, the area in square kilometers of pixels with NDWI above threshold, and valid_fraction, the fraction of AOI pixels that are not NaN. Derive the pixel area from the array’s coordinates, not from a hard-coded number.

In the notebook, use these functions to:

  • Start the cluster as specified in the memory budget section and print the dashboard link.

  • Run search_items for the AOI and period with max_cloud=40, then dedupe_items for tile 10SFH. Print the number of items before and after, and a small table of how many of the original items came from each tile and each processing baseline. In two or three sentences, explain what would have gone wrong in the time series if you had stacked the duplicates.

  • Build the stack with stack_items using only the assets you need (green, nir, scl, that is B03, B08, SCL). Show its repr, its size in GiB (stack.nbytes / 2**30), its chunk size, and its number of chunks, all without computing it.

  • Chain mask_clouds, ndwi, and water_area, then compute the result. This is the only .compute() on imagery in the whole time series step, and the result is two short 1-D arrays.

  • Drop scenes with valid_fraction below 0.9. Plot water area against time as points, with the dropped scenes shown in a different marker, and overlay the monthly median water area of the kept scenes.

  • Pick the dates of the minimum and maximum water area, compute only those two NDWI scenes, and map them side by side with the water mask outlined. Look at the maps critically: if a scene is not believable, find out why from its metadata.

  • Compare your maximum water area to the full-pool area of about 46 km² and explain any difference in two or three sentences.

  • Export the per-scene results to results/water_area_scenes.csv (columns time, water_km2, valid_fraction, kept).

Commit src/utils.py and the current notebook.

Part 2: The agentic precipitation pipeline (25 pts)#

Now use GitHub Copilot in Agent mode, guided by an AGENTS.md you write, to retrieve the precipitation record. The agent adds functions to src/utils.py; your notebook imports them and visualizes the results. Review every change before accepting it. You own this code and these numbers.

2.1 Write your AGENTS.md#

Before prompting the agent, write AGENTS.md. At a minimum it should tell the agent: this is a pixi-only project; all analysis functions go in src/utils.py and must be pure with no plotting; the notebook is for visualization only; do not create new files or edit reflection.md; the domain rules (GHCN-D Parquet on s3://noaa-ghcn-pds read anonymously, read only the years and columns needed, units and quality flags as documented in the GHCN-D readme); the memory rules (the cluster configuration is fixed, functions must return lazy Dask objects unless their result is small by construction, never call .compute() on unreduced data); run tests with pixi run test; and a definition of done. Commit it.

2.2 Functions for the agent to add#

  • load_stations(url): read ghcnd-stations.txt into a pandas DataFrame with at least the station ID, latitude, longitude, and name. (This file is small, so pandas is fine here.)

  • stations_within(stations, lon, lat, radius_km): return the stations within radius_km of a point, with a distance_km column, sorted by distance. Distances must be geodesic (use pyproj.Geod), not computed in degrees.

  • load_ghcn(element, years): return a lazy Dask DataFrame of one element (for example "PRCP") for the given years, read from the Parquet copy with dd.read_parquet and only the columns you need. Inspect the schema (ddf.dtypes) before relying on column names or types.

  • daily_values(ddf, station_ids): filter to the given stations, drop observations that failed quality control, convert values to physical units, and parse dates. Return a lazy Dask DataFrame.

  • completeness(daily, start, end): for each station, the fraction of days between start and end with a valid observation. The result is one number per station, so it is fine to return it computed.

  • monthly_totals(daily_station, min_fraction=0.9): monthly precipitation totals in mm for one station, set to NaN for months in which fewer than min_fraction of days have a valid observation.

2.3 In the notebook, run and visualize#

  • Find all stations within 30 km of the AOI centroid and list them with their distances.

  • Load PRCP for the years needed to cover October 2018 to December 2025 (you need October to December 2018 to complete the first water year, see below), and compute the completeness of every candidate station over that period. Choose the station you will use, and justify the choice in one or two sentences.

  • Filter the daily data to that one station, repartition it to a single partition, and persist it, as in Lecture 14. Then compute monthly totals.

  • Make a two-panel figure with a shared time axis: monthly precipitation (bars) on top and monthly median water area on the bottom.

  • Aggregate both records by water year (October 1st to September 30th, so water year 2023 runs from October 2022 to September 2023): total precipitation, and maximum monthly water area. Scatter one against the other for the water years your imagery fully covers, and label the points by water year.

  • In three to five sentences, interpret the relationship, and name at least two reasons why precipitation at a single station is an imperfect predictor of the area of this reservoir.

  • Export the monthly table to results/folsom_monthly.csv (columns month, water_km2_median, n_scenes, prcp_mm).

Own the data and the memory

These mistakes run without error and produce plausible-looking numbers, or they work on your machine today and crash on someone else’s tomorrow. Watch for the agent: reading the CSVs or a whole year’s archive when only a few columns of one element are needed; calling .compute() on the full DataFrame, or calling it inside a loop instead of once (or once with dask.compute); trusting head() as if it represented the whole dataset (it looks only at the first partition); leaving precipitation in its stored units (check the readme); keeping values that failed quality control; computing distances in degrees; and treating missing days as zero precipitation. Document at least two such mistakes you found and how you corrected them.

Part 3: Memory report (10 pts)#

Add a section to the notebook that measures what your pipeline costs, and save the figure to results/memory.png.

  1. Build, but do not compute, a naive version of the Part 1 stack: all items from a search with no cloud filter and no deduplication; the assets coastal, blue, green, red, rededge1, rededge2, rededge3, nir, nir08, nir09, swir16, swir22, scl (asset keys B01 to B12, B8A, and SCL); resolution=10; and no dtype argument. Call stackstac.stack directly for this; there is no need to convert it to reflectance. Report its number of time steps, its size in GiB, and its number of chunks next to the same numbers for your Part 1 stack, in one small table.

  2. Wrap the Part 1 .compute() and the Part 2 completeness computation in a MemorySampler from distributed.diagnostics, and plot both memory traces:

    from distributed.diagnostics import MemorySampler

ms = MemorySampler()
with ms.sample("water area"):
    result = water.compute()
ms.plot(align=True)

  3. In three to five sentences: which single choice reduced the size of the stack the most, how the peak cluster memory compares with the 4 GiB budget and with the size of the lazy stack, and why the same reduction would still succeed on the naive stack, which is many times larger than the whole cluster’s memory.

Part 4: Tests and quality assurance (10 pts)#

  1. Add pytest to the project and define a test task so the suite runs with pixi run test.

  2. Write tests in tests/test_utils.py on small synthetic inputs where you know the answer. Tests must not touch the network. Include at least: ndwi on a hand-made array with known values; mask_clouds sets a pixel of SCL class 9 to NaN and keeps a pixel of class 4; water_area returns the expected area and valid fraction for a small array with a known number of water pixels at 20 m; dedupe_items keeps the right item from a few hand-built pystac.Item objects that differ in tile and processing baseline; to_reflectance gives the same reflectance for two hand-made scenes of the same surface, one with baseline 03.00 and one with baseline 04.00 (DNs 1000 higher), and leaves scl unchanged; stations_within keeps a point 10 km away and drops one 50 km away; and monthly_totals converts units correctly and returns NaN for a month with too many missing days (build the input with dd.from_pandas).

  3. The agent may help write tests, but verify each test is meaningful: it must fail if the function were wrong. Iterate until pixi run test passes.

Part 5: Reflection (10 pts)#

Write reflection.md yourself, without AI assistance (300 to 500 words):

  • Memory: what did the memory report show you, and if you had to run the same analysis for every reservoir in California over 40 years, which parts of your pipeline would you change first?

  • The agent’s mistakes: describe the incorrect data handling or memory waste the agent produced, how you noticed it, and how writing Part 1 by hand helped you catch it.

  • Retrieval design: GHCN-D is also available partitioned by station (parquet/by_station/). For which step of Part 2 would the by-station layout have been cheaper, and for which step did you need the by-year layout? Also, what does the valid_fraction filter fail to protect you from?

  • Context engineering: how did your AGENTS.md shape the agent’s behavior, and what would you tighten next time?

Deliverables#

Committed to your private geog313-assignments repository under assignment-6/, following the project layout above:

  • pixi.toml, pixi.lock, .gitignore.

  • src/utils.py with all analysis functions, and src/__init__.py.

  • notebooks/analysis.ipynb, executed, with all figures rendered on GitHub, importing from src.utils and containing visualization and measurement only.

  • tests/test_utils.py, with pixi run test passing.

  • results/water_area_scenes.csv, results/folsom_monthly.csv, and results/memory.png.

  • AGENTS.md, README.md, and reflection.md (written without AI).

  • A commit history showing incremental, reviewed development.