Skip to content

kentstephen/crop-data-layer-2021-deckgl-raster

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 

History

22 Commits
 
 
 
 
 
 
 
 
 
 
 
 

Repository files navigation

CDL via deck.gl-raster

This README and the code was generated by Claude wtih my direction. I am including the .claude/memory/ files to add context if someone finds that helpful. Also please check the latest from the Dev Seed repos listed below before you start developing.

Browser-side, viewport-driven, per-class-pickable rendering of USDA's Cropland Data Layer (CDL) from Microsoft Planetary Computer. No tile server, no derived data, no hosting — the browser opens each per-tile COG over HTTP Range and DevSeed's MosaicLayer fans out byte-range reads per visible map tile.

+-------------+   STAC items   +---------------+   signed Range reads
| MPC STAC    |---------------▶| browser app   |◀--+---------------+
| usda-cdl    |                | (deck.gl-     |   | Azure blobs   |
+-------------+                | raster +      |   | (per-tile     |
                               | maplibre)     |   |  COGs +       |
+-------------+   class table  |               |   |  overviews)   |
| GEE STAC    |---------------▶|               |   +---------------+
| catalog     |                +-------+-------+
+-------------+                        │
                                       │ click any pixel
                                       ▼
                                  out-of-band
                                  Range read +
                                  proj4 lon/lat -> CDL pixel
                                  -> class code + crop name

Why this exists

USDA NASS publishes CDL annually as zipped CONUS-wide GeoTIFFs. Microsoft Planetary Computer mirrors that as a STAC collection (usda-cdl) with ~1095 per-tile Cloud-Optimized GeoTIFFs per year, each ~90×90 km in EPSG:5070 Albers with built-in overview pyramids. The data has every property you want for a fast browser viewer (Range-readable, COG-tiled, paletted uint8, embedded colormap) — except it's sharded, not one global COG. Server-side mosaic tile APIs exist but lose per-pixel class codes (you only get the rendered RGB).

DevSeed has the JS side of this solved in deck.gl-raster: MosaicLayer indexes many COGs spatially (Flatbush), fans out viewport-driven byte-range reads per source, and uses each source's own overview pyramid. The naip-mosaic example demonstrates the pattern against MPC's NAIP collection. This repo ports the pattern to CDL with three additions:

  • a 256-entry palette LUT from the GEE STAC catalog (single-band paletted uint8 rendered via a three-module integer pipeline: CreateTextureUintFilterCategoryPaletteColormap)
  • out-of-band per-pixel picking (lon/lat → EPSG:5070 → geotiff.index()fetchTile → read one byte → palette.names[code])
  • a viewport-strict crop dashboard with category filtering and acreage math

Note (May 2026): the renderer was ported from the original r8unorm + normalize-back-to-index hack to the r8uint + usampler2D integer pipeline from the latest examples/land-cover in developmentseed/deck.gl-raster (h/t Kyle Barron). The category filter is now a separate 256-byte boolean LUT texture instead of an alpha-channel mutation on the palette. Same behavior, cleaner shader path, no more snap-to-texel-center math.


Quick start

cd web
npm install

# One-time: build the 256-entry palette LUT from the GEE STAC catalog.
uv run scripts/gen_palette.py

# Re-run when SAS tokens expire (~1 hour). Defaults to 2021 / CONUS.
uv run scripts/gen_stac.py
# uv run scripts/gen_stac.py --year 2020 --bbox -125 24 -66 50

npm run dev   # http://localhost:5454

The two Python scripts are the only Python in the repo. They use PEP-723 inline dependency metadata so uv run handles everything — no separate project install, no venv to manage.

Why Python at all

gen_stac.py pre-bakes SAS-signed hrefs into a static JSON. That sidesteps MPC's sign-endpoint rate limit at CONUS zoom (1095 parallel sign requests trigger 429 cascades). It's one Python step per ~1 hour of token TTL; the rest of the app is pure JS/TS.

gen_palette.py fetches the CDL class table from the Google Earth Engine STAC catalog and bakes it into a JSON the browser imports statically. Only needs to re-run if GEE updates the class list (rare).


What the app does

  • CONUS mosaic at full 30 m resolution. Pan/zoom anywhere; overviews are used at low zoom (the inner COGLayer picks the right pyramid level per source based on screen resolution).
  • Crops-in-viewport dashboard (top-left, collapsible) with color swatch, crop name, acreage, pixel count, and percentage. Recomputes on every pan/zoom and as new tiles land. Acreage from each tile's affine transform: pixels × |a| × |e| × 0.000247105 (m²-to-acres).
  • Category filter with checkboxes for 14 USDA-style groups (Field crops, Cereals, Oilseeds, Legumes, Forage, Fruits/Tree Nuts, Vegetables, Developed, Forest, Shrubland, Pasture, Wetlands, Water, Barren), plus an "Other / fallow / no-data" bucket. Show all / Hide all buttons. Filter affects both the rendered map and the dashboard counts.
  • Click any pixel → bottom-center popup with crop name, CDL class code, and lon/lat.
  • Basemap labels render on top of the raster — the deck.gl layer is inserted before the first symbol layer in whatever maplibre style is loaded (Voyager by default).

Repo layout

.
├── README.md
└── web/
    ├── package.json                vite + react 19 + maplibre + DevSeed packages
    ├── index.html
    ├── vite.config.ts              port 5454; build target esnext (TLA in workers)
    ├── tsconfig.json
    ├── scripts/
    │   ├── gen_stac.py             pystac-client + planetary-computer.sign_inplace
    │   └── gen_palette.py          GEE STAC -> 256-entry RGBA + class names
    └── src/
        ├── main.tsx                React root
        ├── App.tsx                 map, MosaicLayer, dashboard, picking, filter
        ├── proj.ts                 EPSG:5070 + lon/lat <-> Albers transforms
        ├── stats.ts                per-tile histograms + viewport-strict aggregation
        ├── pick.ts                 out-of-band click -> CDL class code + crop name
        ├── cdlShaders.ts           r8uint integer pipeline: CreateTextureUint + FilterCategory + PaletteColormap
        ├── categories.ts           CDL code -> USDA category grouping
        ├── minimal_stac.json       generated, GITIGNORED (SAS-signed, ~1hr TTL)
        └── palette.json            generated (refresh only on GEE updates)

Architecture

Stack

Layer What Notes
Storage Azure blobs (landcoverdata.blob.core.windows.net/usda-cdl/) Per-tile COGs, ~90×90 km in EPSG:5070, embedded overview pyramids + USDA colormap
Catalog MPC STAC (/api/stac/v1) usda-cdl collection; SAS tokens via /api/sas/v1/sign (rate-limited)
Pre-bake pystac-client + planetary-computer.sign_inplace Run hourly to refresh SAS tokens
COG reader @developmentseed/geotiff Range requests, TIFF parser, web-worker decoders (LERC / Deflate / Zstd)
Mosaic @developmentseed/deck.gl-geotiff MosaicLayer (Flatbush spatial index) + COGLayer (pyramid-aware tile fetching)
Render @developmentseed/deck.gl-raster RasterTileLayer, shader-module render pipeline, in-shader EPSG:5070 → Web Mercator
GPU @luma.gl/core + WebGL2/WebGPU r8uint class texture (sampled via usampler2D) + 256×1 RGBA colormap + 256×1 r8 filter LUT
Map maplibre-gl + react-map-gl/maplibre + @deck.gl/mapbox MapboxOverlay interleaves the deck layer beneath basemap labels
Reproject proj4 (TS) + @developmentseed/proj Picking + per-tile bbox projection (TS); GeoTIFF CRS resolution (deck.gl-raster)

Render path (per frame)

  1. MosaicLayer asks: which of the 1095 sources intersect the current viewport bbox? (Flatbush.)
  2. For each: our getSource returns a cached GeoTIFF instance (or opens one — 16 KB header Range read).
  3. renderSource returns a COGLayer per source. The COGLayer's inner RasterTileLayer picks the right overview level for current screen resolution.
  4. For each visible tile: our getTileData does the Range fetch + LERC/Deflate decode, records the tile's histogram + lon/lat bbox in stats.ts, uploads as r8uint to the GPU.
  5. renderTile returns a three-module integer pipeline: [CreateTextureUint, FilterCategory, PaletteColormap]. The first samples the r8uint tile via usampler2D into an ivec4 icolor. The second texelFetches a 256-byte boolean LUT and discards fragments whose class isn't selected. The third texelFetches the 256×1 RGBA colormap and discards alpha=0 (background).
  6. deck.gl composites into maplibre's canvas, beneath the first symbol-layer in the active style.

Picking

Deck.gl's normal picking would only give us the rendered RGB color of a pixel. We do an out-of-band read instead:

  1. maplibre onClick(lng, lat).
  2. Linear bbox-scan over 1095 items → which source COG covers the point.
  3. getCachedGeoTIFF(href) → already-opened GeoTIFF.
  4. proj4(lng, lat) → EPSG:5070 meters.
  5. geotiff.index(x, y) → full-res (row, col).
  6. Derive (tx, ty, localRow, localCol) within the source's internal tile grid.
  7. geotiff.fetchTile(tx, ty) — usually a cache hit if currently visible.
  8. Read one byte → CDL class code → palette.names[code] → popup.

Dashboard

Every tile uploaded to the GPU also feeds a 256-bucket histogram into module-level state (stats.ts), keyed by source and overview level (a switch in overview level for a source drops the prior level's tiles — no double-counting). Each entry stores the tile's lon/lat bbox (4 corners projected from EPSG:5070).

aggregateInViewport(bbox) sums histograms across tiles whose bbox intersects the viewport. Acreage per class = pixel_count × |transform[0]| × |transform[4]| × 0.000247105. Holds at any overview level because pixel area scales with overview.

Granularity is per-tile, not per-pixel: edge tiles that overlap the viewport contribute their full pixel count. Good enough for "what dominates this view"; for survey-grade acreage you'd add per-pixel viewport clipping.

Filter

CDL's 134 classes are grouped into 14 USDA-style categories in categories.ts plus an "Other" bucket for anything unclaimed. When the user toggles a category, we rebuild a 256-byte boolean LUT (255 at every active class code, 0 elsewhere) and upload it as a 256×1 r8unorm texture. The FilterCategory shader module texelFetches the LUT using the integer class code and discards the fragment if the entry is 0. The colormap texture is never touched — palette and filter are decoupled. Categories overlap deliberately (soybeans are both an oilseed and a legume in USDA's classifications) and the UI shows an indeterminate state when that produces partial overlap.


Known footguns

  • SAS tokens expire ~1 hour after gen_stac.py runs. Symptom: 403 / "Server failed to authenticate the request" in the network tab. Fix: re-run the script.
  • MPC's usda-cdl collection caps at 2021 despite USDA publishing through 2024. The MPC mirror is stale.
  • Bilinear filtering on class codes produces garbage. Both textures (class index and palette LUT) use nearest filtering. A "class 1.5" pixel is nonsense.
  • At CONUS zoom (z<5), all 1095 sources are in view. Each contributes its coarsest overview tile. Slow first paint but not broken — the layer's pyramid-aware loading is doing exactly the right thing for the data shape; there's no global pyramid because the data isn't one COG.

What's not (yet) shipped

  • GitHub Pages deployment. Blocked by the 1-hour SAS TTL. Options: runtime sign helper (rate-limit risk; requires a concurrency limiter + retries), or a daily GitHub Action cron (broken between regens). Neither is clean.
  • Per-pixel viewport clipping in stats. Current granularity is per-tile.

Credits

The heavy lifting is all DevSeed:

About

Visualizing USDA Cropland Data Layer from Microsoft Planetary Computer in Lonboard (client-side COG rendering via async-geotiff + obstore)

Resources

Stars

8 stars

Watchers

0 watching

Forks

Releases

No releases published

Packages

 
 
 

Contributors