Workflow for Antarctic Grounding Zone Acceleration

[yueyiche]

Create a jupyter notebook. The goal of this notebook is to use similar tools in issue #25 to look at how the velocity in Antarctic glaciers' grounding zones has been changing with time. Grounding zone here is defined by the coverage of a set of grounding lines (grids on the line plus buffer). We are going to accomplish this with the following workflow.

1. Fetch data through earthaccess. Grounding line data should come from MEaSUREs Antarctic Grounding Line from Differential Satellite Radar Interferometry, Version 2 (DOI: 10.5067/IKBWW4RYHF1Q). Annual ice surface velocity data should come from MEaSUREs Annual Antarctic Ice Velocity Maps, Version 1 (DOI: 10.5067/9T4EPQXTJYW9). If we already have these exact data downloaded, don't download them again, and just use what we have already.
2. Create grounding zone for all Antarctic glaciers for each year. Apply mortie line function to each pieces of grounding line in each year, with a resolution of around 50 m grid. The grids on all the lines plus the buffer of all the lines become our grounding zone. We also create a master grounding zone, where we consider grounding lines from all years, and use mortie on them, so it contains all the grids on all the lines and the buffer for all the lines from 1992 to 2025 (the whole time range for the grounding line data).
3. Intercept annual grounding zone ice velocities. For each year available in the velocity dataset, intercept the grounding zone with the velocity map, and assign the great circle nearest neighbor velocity value to the grounding zone mortie grids. This way we get annual maps of grounding zone velocities for Antarctica. If the grid cell of the certain year did not have an according velocity at the baseline year, map this grid to the 2000-2001 velocity map, and use that as baseline. If the 2000-2001 velocity map has a gap (NA) for that grid cell, then use the first velocity map that has valid data for that grid cell as baseline.
4. Use master grounding zone to intercept velocity maps of 2005 and 2025, similar to what is done in step 3, but use the same master grounding zone for both the 2005 map and the 2025 map. Calculate the velocity anomaly between the two maps.
5. Plot 1: Antarctica grounding zone velocity map slider. Use the 2000-2001 grounding zone velocity as the baseline, make grounding zone velocity anomaly for each grid cell and for each year. Plot a map of Antarctica with these grounding zone locations and the grids with velocity anomaly, red for acceleration, and blue for slow down. Make a slider for these maps, so I can slide through and see how the velocity anomaly has been changing.
6. Plot 2: Antarctica master grounding zone velocity change between 2005 and 2025. Plot the velocity anomaly calculated in step 4.

That should be all. Any questions?

[espg]

Clarifying Q&A — design decisions

Posting before implementation so the choices are reviewable.

Q1 — Render strategy at all-Antarctica scale

A: Compute grounding lines at ~100 m resolution (HEALPix order 16), apply a k=5 morton buffer (~500 m halo), then clip down to the velocity-match order for plotting. Anomaly is rendered at the match order.

Q2 — Per-cell baseline fallback (Plot 1)

A: For each grounding-zone cell, baseline = nearest-neighbour velocity in the 2000–2001 map. If the cell is NA in 2000–2001, walk forward chronologically through the annual maps (2005, 2006, …) and use the first year that has a valid velocity for that cell. Cells that never have valid data → NaN anomaly. anomaly(year, cell) = velocity(year, cell) − baseline(cell).

Q3 — Nearest-neighbour metric

A: No KD-tree. Use morton matching the way the Denman notebook does: morton-index velocity pixels at the match order; intersect with the grounding-zone cells (also at the match order) by morton key.

Q4 — Implementation defaults (multi-select)

A:

• ✓ Reuse examples/insar_gl_cache/ for the NSIDC-0498 geopackage (no re-download).
• ✓ Symmetric anomaly colormap, fixed scale across all slider years (e.g. vmax = global p95 of |anomaly|, vmin = -vmax). RdBu_r so red = acceleration, blue = slow-down.
• ✓ Treat both NaN and exactly-zero velocity values as NA (NSIDC-0720 sometimes uses 0 as a fill).

Q5 — Match order vs actual velocity resolution (raised during planning)

NSIDC-0720 v1 is 1 km, not ~400 m as I'd initially read into the issue. To keep matching honest (one velocity pixel ↔ one morton cell, no rasterisation tricks), match-and-plot at HEALPix order 13 (~800 m) rather than order 14 (~400 m). Grounding lines remain at order 16 (~100 m) → clipped to order 13 by morton-digit truncation.

A: Use order 13 (~800 m) for the match and the plotted resolution.

---

Updates after follow-up comment (2026-04-15)

"yes, the variable names are VX and VY, and I want the velocity magnitude of V."
"For plot 2 actually, instead of intercepting the 2005 and 2025 annual velocity maps with master grounding line, I want to intercept the 1996-2001 and 2020-2022 multi year velocity maps (from MEaSUREs Multi-year Reference Velocity Maps of the Antarctic Ice Sheet, Version 1, DOI: 10.5067/FB851ZIZYX5O), with the master grounding line."

• Velocity variables: VX, VY, computed magnitude V = sqrt(VX² + VY²).
• Plot 2 source dataset changed from NSIDC-0720 (annual, 1 km) to NSIDC-0761 v1 (multi-year reference, 450 m). This is the same dataset used in the Denman notebook — the four epoch files are already cached in examples/insar_vel_cache/ (incl. antarctica_ice_velocity_1995-2001_450m_v01.1.nc and antarctica_ice_velocity_2020-2022_450m_v01.1.nc).
• Note on years: the comment says "1996-2001" but the file in the dataset is 1995-2001 (the actual epoch boundaries are 1995-07-01 → 2001-06-30). Going with the dataset-native epoch.
• Plot 2 match order is now 14 (~400 m) to match the 450 m source, not 13. Plot 1 (NSIDC-0720, 1 km) stays at order 13. So the per-year gz16_year is preserved and re-clipped to 13 or 14 depending on which plot consumes it.

---

Implementation plan

Notebook: examples/morton_grounding_zone_velocity_anomaly.ipynb

Constants

GL_ORDER         = 16       # ~100 m grounding-line cells
BUFFER_K         = 5        # ~500 m halo around each line
PLOT1_MATCH      = 13       # ~800 m, matches 1 km NSIDC-0720 annual velocity
PLOT2_MATCH      = 14       # ~400 m, matches 450 m NSIDC-0761 multi-year velocity
GLACIERS         = 'all'    # full continent, every Glac_Name
PLOT2_EPOCHS     = ('1995-2001', '2020-2022')

Reused helpers (from PR #26)

• download_if_missing(results, local_path, ext) — earthaccess cache shortcut.
• morton_clip_order(m, src, tgt) — digit-truncation between HEALPix orders.
• shapely_to_latlon(geom) — LineString / MultiLineString → lat/lon arrays.
• mortie.linestring_coverage, mortie.morton_buffer, mortie.geo2mort.

Pipeline

1. Fetch data (via download_if_missing):

   • NSIDC-0498 v2 (grounding lines) — reuse cached geopackage.
   • NSIDC-0720 v1 (annual ice velocity, Plot 1) — search short_name='NSIDC-0720', version='1', download all annual files (one 2000–2001 baseline + annual 2005…). Variable count discovered at search time.
   • NSIDC-0761 v1 (multi-year reference, Plot 2) — already cached; only the 1995-2001 and 2020-2022 epoch files are needed.

2. Per-year grounding zone at order 16 (all glaciers):

   • For each Year group in NSIDC-0498:
     • Unpack every feature's LineString/MultiLineString → lat/lon arrays.
     • linestring_coverage(..., order=16) → list of arrays → union → cells16.
     • morton_buffer(cells16, k=5) → ring. gz16_year = union(cells16, ring).
   • GZ16_BY_YEAR[year] = gz16_year (an int64 array of order-16 morton cells).

3. Master grounding zone at order 16:

   • GZ16_MASTER = unique(concatenate([GZ16_BY_YEAR[y] for y in GZ16_BY_YEAR])).

4. Velocity loading helper (parametrised by file path + match order):

   • Open NetCDF with xr.open_dataset(..., engine='h5netcdf').
   • Compute V = sqrt(VX² + VY²); mask NaN and exactly-zero as NA.
   • Reproject pixel centers from EPSG:3031 → EPSG:4326 (single bulk pyproj call).
   • vel_morton = mortie.geo2mort(lats, lons, order=match_order) per pixel.
   • Build a dict[int64, float32] of morton → speed (drop NA pixels).
   • Returns {'speeds': dict, 'span': (start_year, end_year), 'order': match_order}.

5. Plot 1 — load the annual stack at order 13:

   • For each NSIDC-0720 file → VEL_ANNUAL[year] = load_velocity(file, order=13).
   • Pre-clip per-year and master grounding zones from order 16 → 13.

6. Per-cell baseline computation for Plot 1 (over GZ_MASTER_13):

   for cell in GZ_MASTER_13:
       v = VEL_ANNUAL[2000].speeds.get(cell)
       if v is None:
           for year in sorted(VEL_ANNUAL):
               if year == 2000: continue
               v = VEL_ANNUAL[year].speeds.get(cell)
               if v is not None: break
       BASELINE[cell] = v   # may be NaN if never valid
   

   Vectorised via np.isin / np.searchsorted over sorted morton arrays.

7. Per-year anomaly for Plot 1:

   • For each velocity year Y (excluding 2000–2001):
     • For cells in GZ_BY_YEAR_13[Y] where both VEL_ANNUAL[Y] and BASELINE are valid:
       • anomaly[cell] = speed[cell] − BASELINE[cell]
   • ANOMALY_BY_YEAR[Y] = (cells, anomalies) sparse pair of arrays.

8. Master multi-year anomaly for Plot 2 (NSIDC-0720 2005 vs 2025 annual → NSIDC-0761 1995-2001 vs 2020-2022 multi-year):

   • Load the two NSIDC-0761 epoch files at order 14: VEL_REF_BASELINE = load_velocity('1995-2001', order=14), VEL_REF_RECENT = load_velocity('2020-2022', order=14).
   • Pre-clip master grounding zone from order 16 → 14: GZ_MASTER_14.
   • For cells in GZ_MASTER_14 where both reference epochs have valid velocity:
     • master_anomaly[cell] = VEL_REF_RECENT.speeds[cell] − VEL_REF_BASELINE.speeds[cell]

9. Plot 1 — Antarctic slider (south polar stereographic, full continent):

   • Symmetric RdBu_r (red = acceleration, blue = slow-down).
   • Fixed vmax = np.nanpercentile(|all_anomalies|, 95), vmin = -vmax.
   • Render order-13 cells per year. Two render options (will pick the faster one):
     • A. Rasterise to a 1 km EPSG:3031 grid and imshow. One image per slider tick.
     • B. add_geometries(mort2polygon(cells)). Cleaner edges, slower.
   • ipywidgets.SelectionSlider over the years 2005…2025.
   • Faint GZ_MASTER outline as a backdrop for context.

10. Plot 2 — master GZ multi-year anomaly (2005 vs 2025 NSIDC-0720 → 1995-2001 vs 2020-2022 NSIDC-0761):

    • Same colormap and rendering approach as Plot 1, but at order 14 (~400 m, native to the 450 m source).
    • Single static figure. Title makes the two epoch labels explicit.

Memory & speed estimates

• All-Antarctic GLs: ~5–10 M order-16 cells per year (rough), → ~80 MB int64.
• After clip 16 → 13: ~few hundred K cells per year. After clip 16 → 14: ~few × that.
• Master GZ at order 16: ~1–2 M cells; clipped to 13/14: smaller subsets.
• 21-ish NSIDC-0720 annual files × ~few M pixels each → ~100 M morton lookups total. Done lazily per-file, dict per-year stays modest.
• 2 NSIDC-0761 epoch files for Plot 2 — already on disk.
• Slider re-renders should be sub-second since per-year payload is sparse.

Open follow-ups (will resolve while implementing)

• Confirm the actual NSIDC-0720 variable names (VX/VY/V?). The loader will probe lower-case fallbacks like the Denman notebook does. ✓ Confirmed: VX, VY, magnitude V = sqrt(VX² + VY²).
• Confirm there's a 2025 file in v1 (CMR query in step 4 will tell us). ✓ Confirmed: dataset coverage is 2000-2001 baseline + 2005…2025 annual.
• Decide rasterise-vs-polygons (9A vs 9B) once we see the per-year cell counts.

[yueyiche]

yes, the variable names are VX and VY, and I want the velocity magnitude of V.

For plot 2 actually, instead of intercepting the 2005 and 2025 annual velocity maps with master grounding line, I want to intercept the 1996-2001 and 2020-2022 multi year velocity maps (from MEaSUREs Multi-year Reference Velocity Maps of the Antarctic Ice Sheet, Version 1, DOI: 10.5067/FB851ZIZYX5O), with the master grounding line.