From 0b7cf905a6816909ef56d2ba13593695637b050c Mon Sep 17 00:00:00 2001 From: Thomas Teisberg Date: Mon, 14 Sep 2026 14:39:07 -0700 Subject: [PATCH 01/10] Unify deduplication across scripts --- README.md | 6 +++ bedmap_common.py | 75 +++++++++++++++++++++++++++++++++++ plot_bedmap_availability.py | 52 +++--------------------- plot_combined_availability.py | 61 +++------------------------- plot_survey_density.py | 7 +++- 5 files changed, 98 insertions(+), 103 deletions(-) create mode 100644 bedmap_common.py diff --git a/README.md b/README.md index ca413b4..6426c6e 100644 --- a/README.md +++ b/README.md @@ -18,6 +18,12 @@ Or edit `pyproject.toml` to point `xopr` at your local path. ## Entry points +All BedMap consumers share `bedmap_common.py`, which applies the campaign +exclusion list, keeps one version of campaigns present in both BedMap2 and +BedMap3, and parses temporal metadata (a missing or year-9999 end date falls +back to the start date). Change those rules there, not in individual plots. + + ### `plot_combined_availability.py` Stacked bar chart comparing data availability across BedMap and xOPR by year (2001–2023). Categorizes data as open access (xOPR), committed to release, or not released. diff --git a/bedmap_common.py b/bedmap_common.py new file mode 100644 index 0000000..9839b50 --- /dev/null +++ b/bedmap_common.py @@ -0,0 +1,75 @@ +"""Shared BedMap catalog handling. + +Every BedMap consumer (bar charts and the point-data density loader) should +go through here so that campaign exclusions, BM2/BM3 deduplication, and +temporal metadata fixes are applied in exactly one place. +""" + +import re + +import numpy as np +import pandas as pd +from pyproj import Geod +from xopr.bedmap.query import query_bedmap_catalog + +# Campaigns dropped everywhere. CRESIS_2009_Thwaites is a MUSIC swath product +# (~30 cross-track bed picks per radar record) and the Vostok file is scattered +# points; neither can be chained into flight lines. +EXCLUDE_NAMES = { + "RNRF_2008_Vostok-Subglacial-Lake_AIR_BM2", + "CRESIS_2009_Thwaites_AIR_BM3", +} + +_VERSION_SUFFIX = re.compile(r"_BM[123]$") +_GEOD = Geod(ellps="WGS84") + + +def load_bedmap_catalog(collections=("bedmap2", "bedmap3"), exclude=True): + """Return catalog entries with one row per campaign. + + Applies ``EXCLUDE_NAMES``, keeps the newest BedMap version when a campaign + appears in several (BM3 over BM2 over BM1), and parses temporal metadata + into ``ts``/``te``. A missing or sentinel end date (e.g. year 9999) falls + back to the start date. + + Columns: collection, geometry, name, base_name, institution, + temporal_start, temporal_end, ts, te. + """ + cat = query_bedmap_catalog(collections=list(collections)) + props = pd.DataFrame(cat["properties"].tolist(), index=cat.index) + df = pd.concat( + [cat[["collection", "geometry"]], + props[["name", "institution", "temporal_start", "temporal_end"]]], + axis=1, + ) + if exclude: + df = df[~df["name"].isin(EXCLUDE_NAMES)] + df["base_name"] = df["name"].str.replace(_VERSION_SUFFIX, "", regex=True) + # Suffixes sort BM1 < BM2 < BM3, so keep="last" keeps the newest version. + df = df.sort_values("name").drop_duplicates("base_name", keep="last") + + df["ts"] = pd.to_datetime(df["temporal_start"], format="ISO8601", errors="coerce") + te = pd.to_datetime(df["temporal_end"], format="ISO8601", errors="coerce") + bad_end = te.isna() | (te.dt.year > 2100) + df["te"] = te.where(~bad_end, df["ts"]) + return df.reset_index(drop=True) + + +def kept_source_files(collections=("bedmap1", "bedmap2", "bedmap3")): + """Point-data ``source_file`` values to keep after exclusion and dedup. + + ``source_file`` in the point parquet files matches the catalog ``name``. + """ + return set(load_bedmap_catalog(collections)["name"]) + + +def geod_km(geometry): + """Geodesic length in km of a shapely (Multi)LineString in lon/lat.""" + if geometry is None or geometry.is_empty: + return 0.0 + total = 0.0 + for line in (geometry.geoms if hasattr(geometry, "geoms") else [geometry]): + c = np.asarray(line.coords) + if len(c) > 1: + total += _GEOD.inv(c[:-1, 0], c[:-1, 1], c[1:, 0], c[1:, 1])[2].sum() + return total / 1000.0 diff --git a/plot_bedmap_availability.py b/plot_bedmap_availability.py index baef3b9..ed41c77 100644 --- a/plot_bedmap_availability.py +++ b/plot_bedmap_availability.py @@ -8,11 +8,10 @@ from pathlib import Path -import duckdb import matplotlib.pyplot as plt import pandas as pd -from pyproj import Geod -from shapely import wkt + +from bedmap_common import geod_km, load_bedmap_catalog SCRIPT_DIR = Path(__file__).parent OUT_DIR = SCRIPT_DIR / "outputs" @@ -33,12 +32,6 @@ } COUNTRY_ORDER = ["Other", "China", "Russia", "Germany", "UK", "USA"] -# Excluded datasets (scattered points, not reliable flight lines) -EXCLUDE_NAMES = { - "RNRF_2008_Vostok-Subglacial-Lake_AIR_BM2", - "CRESIS_2009_Thwaites_AIR_BM3", -} - def institution_to_country(name): """Map campaign name prefix to country.""" @@ -46,43 +39,10 @@ def institution_to_country(name): return COUNTRY_MAP.get(prefix, "Other") -def line_km(geometry_wkt): - """Calculate geodesic length of a WKT line geometry in km.""" - geom = wkt.loads(geometry_wkt) - geod = Geod(ellps="WGS84") - total = 0.0 - lines = geom.geoms if hasattr(geom, "geoms") else [geom] - for line in lines: - coords = list(line.coords) - for i in range(len(coords) - 1): - _, _, d = geod.inv(coords[i][0], coords[i][1], - coords[i + 1][0], coords[i + 1][1]) - total += d - return total / 1000.0 - - -# Query bedmap2 and bedmap3 catalogs (skip bedmap1, matching reference) -conn = duckdb.connect() -conn.execute("INSTALL spatial; LOAD spatial; SET enable_progress_bar = false;") -urls = [f"https://data.source.coop/englacial/bedmap/bedmap{v}.parquet" for v in [2, 3]] -query = " UNION ALL ".join( - f"SELECT ST_AsText(geometry) as geom_wkt, name, " - f"temporal_start, temporal_end FROM read_parquet('{u}')" for u in urls -) -df = conn.execute(query).fetchdf() -conn.close() - -# Exclude problematic datasets -df = df[~df["name"].isin(EXCLUDE_NAMES)] - -# Deduplicate: keep BM3 over BM2 for campaigns in both catalogs -df["base_name"] = df["name"].str.replace(r"_BM[123]$", "", regex=True) -df = df.sort_values("name").drop_duplicates(subset="base_name", keep="last") - -# Calculate line-km -df["line_km"] = df["geom_wkt"].apply(line_km) -df["ts"] = pd.to_datetime(df["temporal_start"], format="ISO8601") -df["te"] = pd.to_datetime(df["temporal_end"], format="ISO8601") +# Query bedmap2 and bedmap3 catalogs (skip bedmap1, matching reference). +# Exclusions, BM2/BM3 dedup and date parsing live in bedmap_common. +df = load_bedmap_catalog(["bedmap2", "bedmap3"]) +df["line_km"] = df["geometry"].apply(geod_km) # Distribute every campaign evenly across its year range (matching reference) rows = [] diff --git a/plot_combined_availability.py b/plot_combined_availability.py index d4daa14..af6145d 100644 --- a/plot_combined_availability.py +++ b/plot_combined_availability.py @@ -9,13 +9,12 @@ import argparse from pathlib import Path -import duckdb import matplotlib.pyplot as plt import numpy as np import pandas as pd import xopr -from pyproj import Geod -from shapely import wkt + +from bedmap_common import geod_km, load_bedmap_catalog p = argparse.ArgumentParser(description=__doc__) p.add_argument("--haps", action="store_true", @@ -31,62 +30,12 @@ OUT_DIR = SCRIPT_DIR / "outputs" OUT_DIR.mkdir(exist_ok=True) -# Excluded BedMap datasets (scattered points, not reliable flight lines) -EXCLUDE_NAMES = { - "RNRF_2008_Vostok-Subglacial-Lake_AIR_BM2", - "CRESIS_2009_Thwaites_AIR_BM3", -} - - -def geod_km_wkt(geometry_wkt): - """Geodesic length of a WKT geometry in km.""" - geom = wkt.loads(geometry_wkt) - geod = Geod(ellps="WGS84") - total = 0.0 - lines = geom.geoms if hasattr(geom, "geoms") else [geom] - for line in lines: - coords = list(line.coords) - for i in range(len(coords) - 1): - _, _, d = geod.inv(coords[i][0], coords[i][1], - coords[i + 1][0], coords[i + 1][1]) - total += d - return total / 1000.0 - - -def geod_km(geometry): - """Geodesic length of a shapely geometry in km.""" - if geometry is None or geometry.is_empty: - return 0.0 - geod = Geod(ellps="WGS84") - coords = list(geometry.coords) - total = 0.0 - for i in range(len(coords) - 1): - _, _, d = geod.inv(coords[i][0], coords[i][1], - coords[i + 1][0], coords[i + 1][1]) - total += d - return total / 1000.0 - - # --- BedMap line-km per year (2000-2020) --- (Antarctic only, skip for Greenland) if not args.greenland: print("Querying BedMap catalogs...") - conn = duckdb.connect() - conn.execute("INSTALL spatial; LOAD spatial; SET enable_progress_bar = false;") - urls = [f"https://data.source.coop/englacial/bedmap/bedmap{v}.parquet" for v in [2, 3]] - query = " UNION ALL ".join( - f"SELECT ST_AsText(geometry) as geom_wkt, name, " - f"temporal_start, temporal_end FROM read_parquet('{u}')" for u in urls - ) - bm = conn.execute(query).fetchdf() - conn.close() - - bm = bm[~bm["name"].isin(EXCLUDE_NAMES)] - bm["base_name"] = bm["name"].str.replace(r"_BM[123]$", "", regex=True) - bm = bm.sort_values("name").drop_duplicates(subset="base_name", keep="last") - bm["line_km"] = bm["geom_wkt"].apply(geod_km_wkt) - bm["ts"] = pd.to_datetime(bm["temporal_start"], format="ISO8601") - bm["te"] = pd.to_datetime(bm["temporal_end"], format="ISO8601", errors="coerce") - bm["te"] = bm["te"].fillna(bm["ts"]) + # Exclusions, BM2/BM3 dedup and date parsing live in bedmap_common. + bm = load_bedmap_catalog(["bedmap2", "bedmap3"]) + bm["line_km"] = bm["geometry"].apply(geod_km) bm_rows = [] for _, r in bm.iterrows(): diff --git a/plot_survey_density.py b/plot_survey_density.py index 0348f21..b3a90bb 100644 --- a/plot_survey_density.py +++ b/plot_survey_density.py @@ -88,7 +88,12 @@ def load_bedmap(epsg, max_point_spacing_m=1000, local_cache=True): df = query_bedmap(collections=["bedmap1", "bedmap2", "bedmap3"], columns=["lon", "lat", "source_file", "row"], local_cache=local_cache, show_progress=True) - print(f" {len(df)} BedMap points from {df['source_file'].nunique()} files") + # Drop excluded campaigns and BM2/BM3 duplicates (see bedmap_common) + from bedmap_common import kept_source_files + n_files = df["source_file"].nunique() + df = df[df["source_file"].isin(kept_source_files(["bedmap1", "bedmap2", "bedmap3"]))] + print(f" {len(df)} BedMap points from {df['source_file'].nunique()} files " + f"({n_files - df['source_file'].nunique()} excluded or duplicate files dropped)") df = df.sort_values(["source_file", "row"]) xs, ys = tf.transform(df["lon"].values, df["lat"].values) # Mask transitions between files so we don't connect unrelated points From d3e49be4d4d899a9576170aa3dab75f689683e85 Mon Sep 17 00:00:00 2001 From: Thomas Teisberg Date: Mon, 14 Sep 2026 17:12:25 -0700 Subject: [PATCH 02/10] Add AWI, KOPRI helicopter and xOPR-track sources to all figures - extra_sources.py: AWI ice-sounder tracks from the AWI WFS layer (dated per profile, product duplicates removed) replace the AWI BedMap files; KOPRI/UTIG helicopter surveys (QIceRadar ASE2-ASE6) with confirmed seasons. - opr_tracks.py: real flight geometry from CReSIS per-segment CSVs for collections in bedmap_common.OPR_SUBSTITUTIONS (2009 Twin Otter swath case, xOPR-only 2022-2024 seasons); --verify compares CSV/STAC/BedMap lengths. - bedmap_common: drop BM2 files superseded by renamed BM3 files, substituted files and AWI files; campaign_years() puts dated sources in their season. - Bar charts extended to 2025; combined chart reconciles xOPR only against the CReSIS/NASA/UTIG pool so other providers add on top. - Density and coastal scripts include the extra tracks (--no-extra to disable). Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_019UvrC4evHHkT7qWWL96sD5 --- README.md | 34 ++++- bedmap_common.py | 58 +++++++- .../20260914-add-awi-and-kopri-sources.md | 140 ++++++++++++++++++ coastal_survey_gaps.py | 6 +- extra_sources.py | 127 ++++++++++++++++ opr_tracks.py | 130 ++++++++++++++++ plot_bedmap_availability.py | 24 +-- plot_combined_availability.py | 53 ++++--- plot_survey_density.py | 21 ++- 9 files changed, 545 insertions(+), 48 deletions(-) create mode 100644 claude_plans/20260914-add-awi-and-kopri-sources.md create mode 100644 extra_sources.py create mode 100644 opr_tracks.py diff --git a/README.md b/README.md index 6426c6e..37b8348 100644 --- a/README.md +++ b/README.md @@ -19,14 +19,40 @@ Or edit `pyproject.toml` to point `xopr` at your local path. ## Entry points All BedMap consumers share `bedmap_common.py`, which applies the campaign -exclusion list, keeps one version of campaigns present in both BedMap2 and -BedMap3, and parses temporal metadata (a missing or year-9999 end date falls -back to the start date). Change those rules there, not in individual plots. +exclusion list, drops BM2 files superseded by renamed BM3 files, keeps one +version of campaigns present in both BedMap2 and BedMap3, and parses temporal +metadata (a missing or year-9999 end date falls back to the start date). +Change those rules there, not in individual plots. + +Direct provider sources live in `extra_sources.py` and are added to every +BedMap-based figure (bar charts, survey density, coastal gaps; `--no-extra` +disables them in the map scripts): + +- **AWI** ice-sounder tracks (EMR, UWB, UWBM) from the AWI map server WFS + layer `radar:tracks_3031`, dated per profile, 1994/95 to date. These + replace the AWI BedMap files (`EXTERNAL_SOURCE_PREFIXES`). The layer lists + the 60/600 ns EMR pulses and UWB standard/quicklook products as separate + features; `load_awi_campaigns` keeps one per flight. +- **KOPRI/UTIG helicopter radar** (Amundsen Sea, layers ASE2-ASE6) from the + QIceRadar index (Lindzey 2026, Zenodo 21964546 v0.3.0). Seasons are set in + `KOPRI_SEASONS`; the raw data are not released. +- **xOPR substitutes** (`opr_tracks.py`): collections listed in + `bedmap_common.OPR_SUBSTITUTIONS` are loaded from the CReSIS per-segment + CSV flight tracks (real geometry, not the simplified STAC lines) in place of + BedMap files that cannot be chained into flight lines (the 2009 Thwaites + swath product), or for xOPR-only seasons not in BedMap3 (UTIG COLDEX + 2022/23 and 2023/24, GHOST2 2024/25). The institution comes from the STAC + ``opr:provider`` field. `uv run python opr_tracks.py --verify` compares CSV track, + STAC and BedMap lengths for the mapped collections. + +Downloads are cached under `radar_cache/` (gitignored). Dated sources are +placed in their July-June season year; BedMap campaigns keep the calendar-year +spread of their metadata range (`campaign_years`). ### `plot_combined_availability.py` -Stacked bar chart comparing data availability across BedMap and xOPR by year (2001–2023). Categorizes data as open access (xOPR), committed to release, or not released. +Stacked bar chart comparing data availability across BedMap and xOPR by year (2001–2025). Categorizes data as open access (xOPR), committed to release (AWI, including data released outside xOPR; UTIG 2008+), or not released. ```bash uv run python plot_combined_availability.py diff --git a/bedmap_common.py b/bedmap_common.py index 9839b50..1f142df 100644 --- a/bedmap_common.py +++ b/bedmap_common.py @@ -20,6 +20,33 @@ "CRESIS_2009_Thwaites_AIR_BM3", } +# BM2 files whose points all lie within 50 m of a differently named BM3 file +# (renamed on resubmission; checked at point level 2026-09-14), so the suffix +# dedup below cannot catch them. UTIG_2008_ICECAP is 88 % inside UTIG_2010_ICECAP. +SUPERSEDED_NAMES = { + "RNRF_2003_48RAEap5_AIR_BM2", "RNRF_2004_49RAEap5_AIR_BM2", + "RNRF_2005_50RAEap5_AIR_BM2", "RNRF_2006_51RAEap5_AIR_BM2", + "RNRF_2007_52RAEap5_AIR_BM2", "RNRF_2008_53RAEap5_AIR_BM2", + "RNRF_2007_Mirny-Vostok_AIR_BM2", "RNRF_2006_KV1-area_AIR_BM2", + "UTIG_2008_ICECAP_AIR_BM2", +} + +# xOPR collections loaded from real flight tracks (opr_tracks.py): keys are +# collections, values the BedMap files they replace (points that cannot be +# chained into flight lines) or [] for xOPR-only seasons (UTIG COLDEX 2022/23, +# 2023/24 and GHOST2 2024/25 are not in BedMap3). +OPR_SUBSTITUTIONS = { + "2009_Antarctica_TO": ["CRESIS_2009_Thwaites_AIR_BM3", + "CRESIS_2009_AntarcticaTO_AIR_BM3"], + "2022_Antarctica_BaslerMKB": [], + "2023_Antarctica_BaslerMKB": [], + "2024_Antarctica_GroundGHOST2": [], +} + +# Institutions taken entirely from a direct, dated source (extra_sources.py) +# instead of their BedMap submissions. +EXTERNAL_SOURCE_PREFIXES = ("AWI_",) + _VERSION_SUFFIX = re.compile(r"_BM[123]$") _GEOD = Geod(ellps="WGS84") @@ -27,10 +54,11 @@ def load_bedmap_catalog(collections=("bedmap2", "bedmap3"), exclude=True): """Return catalog entries with one row per campaign. - Applies ``EXCLUDE_NAMES``, keeps the newest BedMap version when a campaign - appears in several (BM3 over BM2 over BM1), and parses temporal metadata - into ``ts``/``te``. A missing or sentinel end date (e.g. year 9999) falls - back to the start date. + Applies ``EXCLUDE_NAMES``, ``SUPERSEDED_NAMES``, ``OPR_SUBSTITUTIONS`` and + ``EXTERNAL_SOURCE_PREFIXES``, keeps the newest BedMap version when a + campaign appears in several (BM3 over BM2 over BM1), and parses temporal + metadata into ``ts``/``te``. A missing or sentinel end date (e.g. year + 9999) falls back to the start date. Columns: collection, geometry, name, base_name, institution, temporal_start, temporal_end, ts, te. @@ -43,7 +71,8 @@ def load_bedmap_catalog(collections=("bedmap2", "bedmap3"), exclude=True): axis=1, ) if exclude: - df = df[~df["name"].isin(EXCLUDE_NAMES)] + df = df[~df["name"].isin(dropped_names()) + & ~df["name"].str.startswith(EXTERNAL_SOURCE_PREFIXES)] df["base_name"] = df["name"].str.replace(_VERSION_SUFFIX, "", regex=True) # Suffixes sort BM1 < BM2 < BM3, so keep="last" keeps the newest version. df = df.sort_values("name").drop_duplicates("base_name", keep="last") @@ -55,6 +84,25 @@ def load_bedmap_catalog(collections=("bedmap2", "bedmap3"), exclude=True): return df.reset_index(drop=True) +def dropped_names(): + """BedMap file names removed everywhere (excluded, superseded, substituted).""" + subst = {n for names in OPR_SUBSTITUTIONS.values() for n in names} + return EXCLUDE_NAMES | SUPERSEDED_NAMES | subst + + +def campaign_years(row): + """Years a campaign row contributes to. + + Rows from dated sources carry a ``season`` (July-June start year) and go + entirely into that year; BedMap rows are spread over the calendar years + of their ``ts``..``te`` metadata range (existing convention). + """ + season = row.get("season") if hasattr(row, "get") else None + if season is not None and not pd.isna(season): + return [int(season)] + return list(range(row["ts"].year, row["te"].year + 1)) + + def kept_source_files(collections=("bedmap1", "bedmap2", "bedmap3")): """Point-data ``source_file`` values to keep after exclusion and dedup. diff --git a/claude_plans/20260914-add-awi-and-kopri-sources.md b/claude_plans/20260914-add-awi-and-kopri-sources.md new file mode 100644 index 0000000..76f79ec --- /dev/null +++ b/claude_plans/20260914-add-awi-and-kopri-sources.md @@ -0,0 +1,140 @@ +# Plan: add AWI radar-sounder tracks and the UTIG/KOPRI helicopter surveys + +Status: implemented 2026-09-14 on branch `add-direct-sources`; results in `claude_notes/20260914_direct_sources_implementation.md`. +Background: `claude_notes/20260914_qiceradar_assessment.md`, +`claude_notes/20260914_direct_provider_sources.md`. + +## Goal + +Bring two sources that the current BedMap + xOPR pipeline misses into all +figures (bar charts, combined availability, survey density, coastal gaps): + +1. AWI ice sounders (EMR, UWB, UWBM) from the AWI map server, dated per + profile, 1994/95 to 2025/26. Replaces the AWI BedMap files (subset, undated). +2. UTIG/KOPRI helicopter radar over the Amundsen Sea (QIceRadar layers + ASE2 to ASE6), 2017/18 to 2025/26, absent from every other source. + +Do it through one new module so every figure picks the data up the same way, +mirroring how `bedmap_common.py` centralised the BedMap rules. + +## Design + +### 1. `extra_sources.py`: one campaign table for all consumers + +`load_extra_campaigns() -> GeoDataFrame` with the same columns +`load_bedmap_catalog` returns (`name, institution, geometry [lon/lat], +ts, te`) plus `source` ("awi_wfs" | "qiceradar"), `season` (July-June start +year), `line_km`, `access` ("released" | "committed" | "none"), `doi`, +`season_confidence` ("dated" | "high" | "low"). One row per profile/segment +(AWI) or per feature (KOPRI); consumers aggregate. + +Also `extra_track_points(spacing_m=250) -> DataFrame[x, y, source_file]` +that densifies the geometries for the density/coastal grids, so +`plot_survey_density.load_bedmap`-style consumers can concatenate it. + +### 2. AWI ingest + +- Fetch `https://maps.awi.de/services/common/radar/ows`, WFS 2.0 + `GetFeature typeNames=radar:tracks_3031 outputFormat=application/json` + (~245 MB). Cache raw GeoJSON under `radar_cache/awi/` (gitignored) and a + slim parquet next to it; record the fetch date in the parquet metadata. +- Keep `radar_system in {EMR, UWB, UWBM}` (AWI's own ice sounders). Drop + ACCU/SNOW/ASIRAS (accumulation/snow radars) and all non-AWI holdings + (PASIN, MCoRDS, HiCARS, MARFA), which are already covered by BedMap/xOPR. +- Deduplicate the two-product listing (EMR60/EMR600, UWB standard/qlook): + derive a profile key from `prof_id` with the radar/product token removed, + keep the longest feature set per key, then verify per season that + summed length / buffered-union length is within 1.05 (it is 2.0-2.2 + before dedup; see `awi_portal_overlap.py`). If the key approach does not + reach 1.05 for a season, fall back to a buffered-union geometry for that + season. +- `season` from `date_time_start` (month >= 7 → year, else year-1). + `access`: "released" if `doi` is a dataset DOI (not the bibliography + 10.1594/PANGAEA.972094), else "committed" (AWI policy, matches the + current chart's treatment of AWI). +- Expected result after dedup: ~540k km, per season e.g. 2001 44k, 2003 + 44k, 2013 47k, 2022 14k, 2023 13k, 2024 23k, 2025 25k (table in the + provider-sources note). Nothing for 2019-2021. +- Remove BedMap `AWI_*` entries when the AWI source is active (add to a + `SUPERSEDED_BY_EXTRA` set in `bedmap_common`, applied in + `load_bedmap_catalog`/`kept_source_files`) to avoid double counting. + Note this changes historical AWI bars, roughly +60k km over 2010-2014 + (flights never submitted to BedMap). + +### 3. UTIG/KOPRI helicopter ingest + +- Download the QIceRadar gpkg from the pinned DOI 10.5281/zenodo.21964546 + (v0.3.0) into `radar_cache/qiceradar/`; read layers ASE2..ASE6 only. + Reproject 3031 → 4326. +- Season table (manual, in the module). KOPRI Araon cruises are annual and + visit Thwaites every other year, most recently 2025/26 (user, 2026-09-14), + which fixes the alternating sequence: + KRT1 2016, ASE2 2017 (Dodson/Getz), KRT2 2018, ASE3 2019 (ANA10B, files + dated 31 Jan-11 Feb 2020, TDR doi:10.18738/T8/730HSL), ASE4 2021 (ANA12, + files dated 31 Jan-4 Feb 2022; the TDR title/date fields say 2023 but the + file names, KML label, Level 2 title and build notes say 2022; TDR + doi:10.18738/T8/3RSCPO), ASE5 2023, ASE6 2025. All `access` = "none": + TDR hosts Level 1B focused radargrams and Level 2 tables, not raw data, + and no release commitment exists. `institution` = "KOPRI", country "Other". +- ~14.8k km total; ~9.7k km within 150 km of Thwaites. +- No overlap with BedMap/xOPR, so nothing to remove. + +### 4. Consumers + +- `plot_bedmap_availability.py` and `plot_combined_availability.py`: + concatenate `load_bedmap_catalog(...)` with `load_extra_campaigns()`. + Extra rows carry a single `season`; put their line-km in that year (no + calendar spread, because they have real dates). Combined chart: AWI + "released" and KOPRI "released" go to the "commitment" bucket unless a + fourth category "released outside xOPR" is added (decision below). + Extend the x-axis to 2025 (decision below). +- `plot_survey_density.py` / `coastal_survey_gaps.py`: when `--source` + is a BedMap variant, append `extra_track_points()` to the point table + before `bin_line_km`; add `--no-extra` to disable. The 1 km gap filter + is irrelevant at 250 m spacing. +- `README.md`: document the new module, the cache, the pinned versions. + +### 5. Verification + +- Per-season AWI km after dedup matches the table in the provider note + within a few percent; buffered-union ratio <= 1.05 per season. +- BedMap AWI files no longer appear in `kept_source_files`; total AWI km + in the bar chart = extra-source AWI km. +- Cell-overlap check (reuse `qiceradar_spatial_overlap.py` machinery): + extra-source km inside xOPR cells stays < 2 %. +- Re-run all four figures; compare summary totals before/after and record + them in `claude_notes/`. + +## Decisions (user, 2026-09-14) + +1. Replace all AWI BedMap entries with the portal source. +2. Extend chart windows to 2025. +3. Three categories. All publicly available AWI data counts as "commitment + to release" (it is not in xOPR); KOPRI helicopter data is "no access". +4. New sources by season year; existing BedMap handling unchanged. +5. Fix the renamed BM2/BM3 duplicates in the same change. +6. Added scope: tooling to substitute xOPR collections for unreliable + BedMap files (the 2009 CReSIS swath case), loading real flight geometry + rather than the simplified STAC geometry, with a length check against + BedMap for every substituted file (expected to agree except 2009). +7. PRIC is out of scope (season labels untrusted). + +## Decisions needed before implementing (resolved above) + +1. Replace all AWI BedMap entries with the WFS source (recommended: dated, + complete, superset) or only add seasons >= 2020 (keeps historical bars + unchanged)? +2. Extend chart windows to 2025 (data exists for 2022/23-2025/26)? +3. Combined chart: keep three categories (AWI/KOPRI releases counted as + "commitment") or add "released outside xOPR"? +4. Season-year assignment for the new sources (recommended) while BedMap + keeps the calendar spread the user chose to leave alone; accept the + mixed convention or revisit 3a of the repo review at the same time. +5. Fix the renamed BM2/BM3 duplicates (section 5 of the QIceRadar note, + ~55k km) in the same `bedmap_common` change? + +## Out of scope + +PRIC/CHINARE lines (season labels untrusted, see the QIceRadar note); +other QIceRadar-only layers (unreleased UTIG SOAR surveys, BAS 1969-1988, +SPRI), which are pre-window or undated. diff --git a/coastal_survey_gaps.py b/coastal_survey_gaps.py index 8defb02..47e3539 100644 --- a/coastal_survey_gaps.py +++ b/coastal_survey_gaps.py @@ -447,6 +447,8 @@ def main(): ) p.add_argument('--source', choices=['bedmap', 'bedmap_local', 'xopr'], default='bedmap_local') + p.add_argument('--no-extra', action='store_true', + help='Do not add extra_sources tracks to the BedMap sources') p.add_argument('--grid-km', type=float, default=30) p.add_argument('--coast-dist-km', type=float, default=20, help='Max distance from coastline [km]') @@ -475,9 +477,9 @@ def main(): # Load survey data print(f'Loading {args.source} data...') if args.source == 'bedmap': - x1, y1, x2, y2 = load_bedmap(reg['epsg'], local_cache=False) + x1, y1, x2, y2 = load_bedmap(reg['epsg'], local_cache=False, include_extra=not args.no_extra) elif args.source == 'bedmap_local': - x1, y1, x2, y2 = load_bedmap(reg['epsg'], local_cache=True) + x1, y1, x2, y2 = load_bedmap(reg['epsg'], local_cache=True, include_extra=not args.no_extra) else: geoms = load_xopr(region_filter='Antarctica') x1, y1, x2, y2 = extract_segments(geoms, reg['epsg']) diff --git a/extra_sources.py b/extra_sources.py new file mode 100644 index 0000000..ee89a60 --- /dev/null +++ b/extra_sources.py @@ -0,0 +1,127 @@ +"""Campaign tracks from direct provider sources, for all figures. + +* AWI ice sounders (EMR, UWB, UWBM) from the AWI map server WFS layer + ``radar:tracks_3031``, dated per profile, 1994/95 to date. Replaces the + AWI BedMap files (see ``bedmap_common.EXTERNAL_SOURCE_PREFIXES``). +* UTIG/KOPRI helicopter radar over the Amundsen Sea from the QIceRadar index + (Lindzey 2026, v0.3.0), layers ASE2-ASE6. Absent from BedMap and xOPR. +* xOPR collections substituted for unusable BedMap files (``opr_tracks``). + +``load_extra_campaigns()`` returns rows shaped like ``load_bedmap_catalog`` +plus ``source``, ``season`` (July-June start year), ``access`` and +``line_km``; ``extra_segments(epsg)`` returns segment endpoints for gridding. +Raw downloads are cached under ``radar_cache/`` (gitignored). +""" + +import urllib.request +from pathlib import Path + +import geopandas as gpd +import numpy as np +import pandas as pd +from pyproj import Transformer + +from bedmap_common import geod_km + +CACHE = Path(__file__).parent / "radar_cache" +AWI_WFS = ("https://maps.awi.de/services/common/radar/ows?service=WFS&version=2.0.0" + "&request=GetFeature&typeNames=radar:tracks_3031&outputFormat=application/json") +AWI_SOUNDERS = {"EMR", "UWB", "UWBM"} # ACCU/SNOW/ASIRAS are accumulation radars +AWI_GENERIC_DOI = "10.1594/PANGAEA.972094" # collection bibliography, not a dataset release +QICERADAR_GPKG = "https://zenodo.org/api/records/21964546/files/qiceradar_antarctic_index.gpkg/content" +# KOPRI Araon cruises are annual and visit Thwaites every other year (last +# 2025/26); ASE3/ASE4 dates confirmed from the TDR file names. +KOPRI_SEASONS = {"ASE2": 2017, "ASE3": 2019, "ASE4": 2021, "ASE5": 2023, "ASE6": 2025} + + +def _fetch(url, path): + if not path.exists(): + path.parent.mkdir(parents=True, exist_ok=True) + urllib.request.urlretrieve(url, path) + return path + + +def _season(ts): + return ts.dt.year.where(ts.dt.month >= 7, ts.dt.year - 1) + + +def load_awi_tracks(): + """All features of the AWI WFS layer (EPSG:3031) with ``km`` added; cached as parquet.""" + pq = CACHE / "awi" / "tracks_3031.parquet" + if pq.exists(): + return gpd.read_parquet(pq) + g = gpd.read_file(_fetch(AWI_WFS, CACHE / "awi" / "tracks_3031.geojson")).set_crs(3031, allow_override=True) + g = g[g.geom_type.isin(["LineString", "MultiLineString"])].copy() + g["km"] = g.to_crs(4326)["geometry"].apply(geod_km) + g.to_parquet(pq) + return g + + +def load_awi_campaigns(): + """AWI ice-sounder tracks, one row per profile segment, duplicates removed.""" + g = load_awi_tracks() + g = g[g["radar_system"].isin(AWI_SOUNDERS)].copy() + # UWB/UWBM list a 'standard' and a 'qlook' product per granule: keep one. + qlook = g["quick_looks"].fillna("").str.contains("_qlook") + key = g["prof_id"] + "|" + g["date_time_start"].astype(str) + g = g[~(qlook & key.isin(key[~qlook]))] + # EMR lists the 60 ns and 600 ns pulse as separate profiles of the same + # flight (prof_id 'antrYYYY_..._EMR60[0]_...' or an 8-digit flight number + # whose third digit is the product): keep the longer product per flight. + flight = (g["prof_id"].str.replace(r"_EMR600?_", "_EMR_", regex=True) + .str.replace(r"^(\d{4})\d(\d{3})$", r"\1_\2", regex=True)) + best = g.groupby([flight, g["radar_name"]])["km"].sum().groupby(level=0).idxmax().str[1] + g = g[g["radar_name"].values == flight.map(best).values] + ts = pd.to_datetime(g["date_time_start"], utc=True) + out = gpd.GeoDataFrame({ + "name": "AWI_" + g["season"].str.replace(" ", "_"), + "institution": "AWI", "source": "awi_wfs", + "season": _season(ts).astype(int), "ts": ts, "te": pd.to_datetime(g["date_time_end"], utc=True), + "radar_system": g["radar_system"], "doi": g["doi"], + "released": g["doi"].notna() & ~g["doi"].str.contains(AWI_GENERIC_DOI, na=False), + "access": "committed", # AWI policy; released data are not in xOPR + "line_km": g["km"], + }, geometry=g.to_crs(4326)["geometry"].values, crs="EPSG:4326") + return out.reset_index(drop=True) + + +def load_kopri_campaigns(): + """UTIG/KOPRI helicopter radar tracks (QIceRadar ASE layers).""" + gpkg = _fetch(QICERADAR_GPKG, CACHE / "qiceradar" / "qiceradar_antarctic_index.gpkg") + parts = [] + for layer, season in KOPRI_SEASONS.items(): + g = gpd.read_file(gpkg, layer=layer).to_crs(4326) + g = g[g.geom_type == "LineString"] + parts.append(gpd.GeoDataFrame({ + "name": f"KOPRI_{layer}", "institution": "KOPRI", "source": "qiceradar", + "season": season, "access": "none", "line_km": g["geometry"].apply(geod_km), + }, geometry=g["geometry"].values, crs="EPSG:4326")) + return pd.concat(parts, ignore_index=True) + + +def load_extra_campaigns(include_opr=True): + parts = [load_awi_campaigns(), load_kopri_campaigns()] + if include_opr: + from opr_tracks import opr_campaigns + parts.append(opr_campaigns()) + return pd.concat(parts, ignore_index=True) + + +def extra_segments(epsg, campaigns=None): + """Segment endpoints (x1, y1, x2, y2) in ``epsg`` for gridding line-km.""" + g = load_extra_campaigns() if campaigns is None else campaigns + tf = Transformer.from_crs("EPSG:4326", epsg, always_xy=True) + x1, y1, x2, y2 = [], [], [], [] + for geom in g["geometry"]: + for line in (geom.geoms if hasattr(geom, "geoms") else [geom]): + c = np.asarray(line.coords) + if len(c) < 2: + continue + x, y = tf.transform(c[:, 0], c[:, 1]) + x1.append(x[:-1]); y1.append(y[:-1]); x2.append(x[1:]); y2.append(y[1:]) + return tuple(np.concatenate(a) for a in (x1, y1, x2, y2)) + + +if __name__ == "__main__": + c = load_extra_campaigns() + print(c.groupby(["source", "institution", "season"])["line_km"].sum().round(0).to_string()) diff --git a/opr_tracks.py b/opr_tracks.py new file mode 100644 index 0000000..53493d5 --- /dev/null +++ b/opr_tracks.py @@ -0,0 +1,130 @@ +"""Real flight tracks for xOPR collections, from CReSIS per-segment CSVs. + +The STAC catalog geometry is simplified, so for line-km and gridding we use +``https://data.cresis.ku.edu/data/rds//csv/Data_.csv`` +(full-rate LAT/LON/TIME/FRAME per record, ~3 MB per segment), cached under +``radar_cache/opr_csv/``. Used for the collections in +``bedmap_common.OPR_SUBSTITUTIONS`` (BedMap files that cannot be chained into +lines, and xOPR-only seasons); ``python opr_tracks.py --verify`` compares +track lengths with the BedMap files they replace or duplicate. +""" + +import re +import urllib.request +from pathlib import Path + +import geopandas as gpd +import numpy as np +import pandas as pd +from pyproj import Geod +from shapely.geometry import LineString, MultiLineString + +from bedmap_common import OPR_SUBSTITUTIONS, geod_km + +CSV_URL = "https://data.cresis.ku.edu/data/rds/{collection}/csv/" +CACHE = Path(__file__).parent / "radar_cache" / "opr_csv" +GAP_M = 2000 # split a segment where consecutive records are further apart +_GEOD = Geod(ellps="WGS84") + + +def _fetch(url, path): + if not path.exists(): + path.parent.mkdir(parents=True, exist_ok=True) + urllib.request.urlretrieve(url, path) + return path + + +def list_segment_csvs(collection): + """Segment CSV file names on the CReSIS server (cached listing).""" + idx = _fetch(CSV_URL.format(collection=collection), CACHE / collection / "index.html") + return sorted(set(re.findall(r'href="(Data_\d{8}_\d{2}\.csv)"', idx.read_text()))) + + +def _segment_line(df): + """Chain records into a (Multi)LineString, breaking at gaps > GAP_M.""" + lon, lat = df["LON"].to_numpy(), df["LAT"].to_numpy() + d = _GEOD.inv(lon[:-1], lat[:-1], lon[1:], lat[1:])[2] + breaks = np.flatnonzero(d > GAP_M) + 1 + parts = [LineString(c) for c in np.split(np.column_stack([lon, lat]), breaks) if len(c) > 1] + return MultiLineString(parts) if len(parts) > 1 else parts[0] + + +def load_opr_tracks(collection, only_stac_segments=True): + """One row per segment: collection, segment, date, season, provider (STAC + ``opr:provider``), geometry (lon/lat), line_km.""" + stac_segments, provider = None, "CRESIS" + if only_stac_segments: + import xopr + items = xopr.OPRConnection().query_frames(collections=[collection], exclude_geometry=True) + stac_segments = {f"{p['opr:date']}_{int(p['opr:segment']):02d}" for p in items["properties"]} + provider = pd.Series([p.get("opr:provider", "cresis") for p in items["properties"]]).mode()[0].upper() + rows = [] + for name in list_segment_csvs(collection): + seg = name[5:-4] + if stac_segments is not None and seg not in stac_segments: + continue + df = pd.read_csv(_fetch(CSV_URL.format(collection=collection) + name, CACHE / collection / name), + usecols=["LAT", "LON"]) + df = df[(df["LAT"].abs() <= 90) & (df["LON"].abs() <= 180)] + if len(df) < 2: + continue + date = pd.Timestamp(seg[:8]) + rows.append({"collection": collection, "segment": seg, "date": date, + "season": date.year if date.month >= 7 else date.year - 1, + "provider": provider, "geometry": _segment_line(df)}) + g = gpd.GeoDataFrame(rows, geometry="geometry", crs="EPSG:4326") + g["line_km"] = g["geometry"].apply(geod_km) + return g + + +def opr_campaigns(): + """Campaign rows (one per collection) for the substituted xOPR collections.""" + rows = [] + for coll in OPR_SUBSTITUTIONS: + t = load_opr_tracks(coll) + season = int(t["season"].mode().iloc[0]) + inst = t["provider"].iloc[0] + rows.append({"name": f"{inst}_{coll}", "institution": inst, "source": "opr", + "season": season, "access": "open", "line_km": t["line_km"].sum(), + "geometry": MultiLineString([p for g in t["geometry"] + for p in (g.geoms if hasattr(g, "geoms") else [g])])}) + return gpd.GeoDataFrame(rows, geometry="geometry", crs="EPSG:4326") + + +# BedMap files that correspond to xOPR collections (from QIceRadar's +# bedmap_labels.py); used only by --verify. +VERIFY_MAP = { + "2004_Antarctica_P3chile": ["NASA_2004_ICEBRIDGE_AIR_BM2"], + "2009_Antarctica_TO": OPR_SUBSTITUTIONS["2009_Antarctica_TO"], + "2009_Antarctica_DC8": ["NASA_2009_ICEBRIDGE_AIR_BM2"], + "2010_Antarctica_DC8": ["NASA_2010_ICEBRIDGE_AIR_BM2"], + "2012_Antarctica_DC8": ["NASA_2012_ICEBRIDGE_AIR_BM2"], + "2013_Antarctica_Basler": ["CRESIS_2013_Siple-Coast_AIR_BM3"], + "2013_Antarctica_P3": ["NASA_2013_ICEBRIDGE_AIR_BM3"], + "2019_Antarctica_GV": ["NASA_2019_ICEBRIDGE_AIR_BM3"], +} + + +def verify(): + """Compare CSV track km with the STAC geometry and the matching BedMap files.""" + import xopr + from bedmap_common import load_bedmap_catalog + cat = load_bedmap_catalog(["bedmap2", "bedmap3"], exclude=False).set_index("name") + conn = xopr.OPRConnection() + for coll, bm_names in VERIFY_MAP.items(): + t = load_opr_tracks(coll) + stac = conn.query_frames(collections=[coll], exclude_geometry=False) + stac_km = stac["geometry"].apply(geod_km).sum() + bm_km = sum(geod_km(cat.loc[n, "geometry"]) for n in bm_names if n in cat.index) + print(f"{coll:26s} csv {t['line_km'].sum():9.0f} km | STAC {stac_km:9.0f} km | " + f"BedMap catalog {bm_km:9.0f} km ({', '.join(bm_names)})") + + +if __name__ == "__main__": + import sys + if "--verify" in sys.argv: + verify() + else: + for coll in OPR_SUBSTITUTIONS: + t = load_opr_tracks(coll) + print(f"{coll}: {len(t)} segments, {t['line_km'].sum():.0f} km, seasons {sorted(t['season'].unique())}") diff --git a/plot_bedmap_availability.py b/plot_bedmap_availability.py index ed41c77..1592b4e 100644 --- a/plot_bedmap_availability.py +++ b/plot_bedmap_availability.py @@ -1,9 +1,10 @@ #!/usr/bin/env python3 """Stacked bar chart of BedMap line-km per year, colored by country. -Recreates the reference figure from the BedMap STAC catalog data. -Uses campaign name years from temporal_start/temporal_end metadata, -distributing multi-year campaigns evenly across their span. +Recreates the reference figure from the BedMap STAC catalog data, plus the +direct sources in extra_sources (AWI tracks, KOPRI helicopter surveys, xOPR +substitutes). BedMap campaigns are spread evenly over their +temporal_start..temporal_end calendar years; dated sources go to their season. """ from pathlib import Path @@ -11,7 +12,8 @@ import matplotlib.pyplot as plt import pandas as pd -from bedmap_common import geod_km, load_bedmap_catalog +from bedmap_common import campaign_years, geod_km, load_bedmap_catalog +from extra_sources import load_extra_campaigns SCRIPT_DIR = Path(__file__).parent OUT_DIR = SCRIPT_DIR / "outputs" @@ -43,24 +45,24 @@ def institution_to_country(name): # Exclusions, BM2/BM3 dedup and date parsing live in bedmap_common. df = load_bedmap_catalog(["bedmap2", "bedmap3"]) df["line_km"] = df["geometry"].apply(geod_km) +df = pd.concat([df, load_extra_campaigns()], ignore_index=True) -# Distribute every campaign evenly across its year range (matching reference) +# Distribute every campaign evenly across its years (see campaign_years) rows = [] for _, r in df.iterrows(): - y_start, y_end = r["ts"].year, r["te"].year - n_years = y_end - y_start + 1 + years = campaign_years(r) country = institution_to_country(r["name"]) - for y in range(y_start, y_end + 1): - rows.append({"year": y, "line_km": r["line_km"] / n_years, + for y in years: + rows.append({"year": y, "line_km": r["line_km"] / len(years), "country": country}) result = pd.DataFrame(rows) -result = result[(result["year"] >= 2000) & (result["year"] <= 2020)] +result = result[(result["year"] >= 2000) & (result["year"] <= 2025)] # Pivot: line-km per year per country, ensure all years present pivot = result.pivot_table(index="year", columns="country", values="line_km", aggfunc="sum", fill_value=0) -pivot = pivot.reindex(range(2000, 2021), fill_value=0) +pivot = pivot.reindex(range(2000, 2026), fill_value=0) pivot = pivot.reindex(columns=[c for c in COUNTRY_ORDER if c in pivot.columns]) # Plot diff --git a/plot_combined_availability.py b/plot_combined_availability.py index af6145d..369a1b3 100644 --- a/plot_combined_availability.py +++ b/plot_combined_availability.py @@ -1,9 +1,11 @@ #!/usr/bin/env python3 """Combined figure: BedMap vs xOPR data availability. -For 2000-2019, xOPR data is a subset of BedMap. For 2020+, only xOPR exists. -Three categories: "Open access to raw data" (xOPR), "Commitment to release" -(AWI all years + UTIG 2008+), and "Raw data not released" (remainder). +Totals come from the BedMap catalog plus the direct sources in extra_sources +(AWI tracks, KOPRI helicopter surveys, xOPR substitutes). Three categories: +"Open access to raw data" (xOPR), "Commitment to release" (AWI all years, +including data released outside xOPR, + UTIG 2008+), and "Raw data not +released" (remainder, including the KOPRI helicopter data). """ import argparse @@ -14,7 +16,8 @@ import pandas as pd import xopr -from bedmap_common import geod_km, load_bedmap_catalog +from bedmap_common import campaign_years, geod_km, load_bedmap_catalog +from extra_sources import load_extra_campaigns p = argparse.ArgumentParser(description=__doc__) p.add_argument("--haps", action="store_true", @@ -36,25 +39,31 @@ # Exclusions, BM2/BM3 dedup and date parsing live in bedmap_common. bm = load_bedmap_catalog(["bedmap2", "bedmap3"]) bm["line_km"] = bm["geometry"].apply(geod_km) + bm = pd.concat([bm, load_extra_campaigns()], ignore_index=True) bm_rows = [] for _, r in bm.iterrows(): - y_start, y_end = r["ts"].year, r["te"].year - n_years = y_end - y_start + 1 + years = campaign_years(r) prefix = r["name"].split("_")[0] - committed = prefix == "AWI" or (prefix == "UTIG" and y_start >= 2008) - for y in range(y_start, y_end + 1): - bm_rows.append({"year": y, "line_km": r["line_km"] / n_years, - "committed": committed}) + committed = (r.get("access") == "committed" or prefix == "AWI" + or (prefix == "UTIG" and years[0] >= 2008)) + # Data that OPR may also host (CReSIS/NASA, UTIG, substituted xOPR + # collections) is reconciled against the xOPR total below; other + # providers' data is never in xOPR and is added on top. + in_opr_pool = prefix in ("NASA", "CRESIS", "UTIG") or r.get("source") == "opr" + for y in years: + bm_rows.append({"year": y, "line_km": r["line_km"] / len(years), + "committed": committed, "pool": in_opr_pool}) bm_df = pd.DataFrame(bm_rows) - bedmap_yearly = (bm_df.groupby("year")["line_km"].sum() - .reindex(range(2001, 2024), fill_value=0)) - committed_yearly = (bm_df[bm_df["committed"]].groupby("year")["line_km"].sum() - .reindex(range(2001, 2024), fill_value=0)) + _yr = lambda d: d.groupby("year")["line_km"].sum().reindex(range(2001, 2026), fill_value=0) + pool_yearly = _yr(bm_df[bm_df["pool"]]) + other_yearly = _yr(bm_df[~bm_df["pool"]]) + committed_yearly = _yr(bm_df[bm_df["committed"]]) else: - bedmap_yearly = pd.Series(0, index=range(2001, 2024)) - committed_yearly = pd.Series(0, index=range(2001, 2024)) + pool_yearly = pd.Series(0, index=range(2001, 2026)) + other_yearly = pd.Series(0, index=range(2001, 2026)) + committed_yearly = pd.Series(0, index=range(2001, 2026)) # --- xOPR line-km per year (Antarctic + Greenland) --- print("Querying xOPR catalog...") @@ -82,19 +91,21 @@ def _yearly_sum(rows): if not rows: - return pd.Series(0, index=range(2001, 2024), dtype=float) + return pd.Series(0, index=range(2001, 2026), dtype=float) return (pd.DataFrame(rows).groupby("year")["line_km"].sum() - .reindex(range(2001, 2024), fill_value=0)) + .reindex(range(2001, 2026), fill_value=0)) opr_yearly = _yearly_sum(opr_rows) # Greenland xOPR data is not in BedMap, so add it to the BedMap total greenland_yearly = _yearly_sum(greenland_rows) # --- Combine: open access / commitment / not released --- -# Greenland xOPR is additional to BedMap (Antarctic only), so add it to totals -years = np.arange(2001, 2024) +# Antarctic xOPR data overlaps the CReSIS/NASA/UTIG pool, so take the larger of +# the two; providers outside OPR (AWI, BAS, ...) and Greenland xOPR add on top. +years = np.arange(2001, 2026) open_km = opr_yearly.values -total_km = np.maximum(bedmap_yearly.values + greenland_yearly.values, open_km) +total_km = (np.maximum(pool_yearly.values, open_km - greenland_yearly.values) + + other_yearly.values + greenland_yearly.values) # Committed data not yet in xOPR (cap at remaining BedMap after removing xOPR) remaining = total_km - open_km commit_km = np.minimum(committed_yearly.values, remaining) diff --git a/plot_survey_density.py b/plot_survey_density.py index b3a90bb..fb08ebc 100644 --- a/plot_survey_density.py +++ b/plot_survey_density.py @@ -67,8 +67,8 @@ } -def load_bedmap(epsg, max_point_spacing_m=1000, local_cache=True): - """Load BedMap2+3 point data and return projected segment endpoints. +def load_bedmap(epsg, max_point_spacing_m=1000, local_cache=True, include_extra=True): + """Load BedMap point data (plus extra_sources tracks) as projected segment endpoints. Parameters ---------- @@ -76,6 +76,9 @@ def load_bedmap(epsg, max_point_spacing_m=1000, local_cache=True): Target CRS (e.g. "EPSG:3031"). max_point_spacing_m : float Discard segments longer than this (gap between survey points). + include_extra : bool + Append the direct-source tracks from extra_sources (AWI, KOPRI, xOPR + substitutes), which replace the corresponding BedMap files. Returns ------- @@ -100,7 +103,13 @@ def load_bedmap(epsg, max_point_spacing_m=1000, local_cache=True): same_file = df["source_file"].values[:-1] == df["source_file"].values[1:] dists = np.sqrt(np.diff(xs)**2 + np.diff(ys)**2) ok = same_file & (dists < max_point_spacing_m) & (dists > 0) - return xs[:-1][ok], ys[:-1][ok], xs[1:][ok], ys[1:][ok] + segs = [xs[:-1][ok], ys[:-1][ok], xs[1:][ok], ys[1:][ok]] + if include_extra: + from extra_sources import extra_segments + ex = extra_segments(epsg) + print(f" + {len(ex[0])} segments from extra_sources") + segs = [np.concatenate([a, b]) for a, b in zip(segs, ex)] + return tuple(segs) def load_xopr(region_filter=None): """Load geometries from xOPR STAC catalog. @@ -175,6 +184,8 @@ def bin_line_km(x1, y1, x2, y2, grid_m, max_extent): default="antarctica") p.add_argument("--zoom", choices=list(ZOOM_REGIONS.keys()), help="Zoom to a named sub-region") + p.add_argument("--no-extra", action="store_true", + help="Do not add extra_sources tracks to the BedMap sources") p.add_argument("--grid-km", type=float, default=30) p.add_argument("--target-spacing", type=float, default=None, help="Target survey spacing [km]; uses diverging PiYG colorscale") @@ -191,11 +202,11 @@ def bin_line_km(x1, y1, x2, y2, grid_m, max_extent): if args.source == "bedmap": if args.region == "greenland": print(" Warning: BedMap is Antarctic-only, map will be empty.") - x1, y1, x2, y2 = load_bedmap(reg["epsg"], local_cache=False) + x1, y1, x2, y2 = load_bedmap(reg["epsg"], local_cache=False, include_extra=not args.no_extra) elif args.source == "bedmap_local": if args.region == "greenland": print(" Warning: BedMap is Antarctic-only, map will be empty.") - x1, y1, x2, y2 = load_bedmap(reg["epsg"], local_cache=True) + x1, y1, x2, y2 = load_bedmap(reg["epsg"], local_cache=True, include_extra=not args.no_extra) else: region_filter = {"antarctica": "Antarctica", "greenland": "Greenland"}[args.region] From a65f728d157d39231d0cd1658405cdf7814f88a1 Mon Sep 17 00:00:00 2001 From: Thomas Teisberg Date: Mon, 14 Sep 2026 19:12:59 -0700 Subject: [PATCH 03/10] Bar charts use gap-filtered point line-km; shared 1.5 km gap rule - bedmap_common: load_bedmap_points/point_segments hold the one gap rule (MAX_POINT_SPACING_M, raised 1 km -> 1.5 km to keep near-1 km airborne surveys while still dropping traverses and 25 km-spaced GANOVEX files); campaign_point_km() caches per-file km in bedmap_campaign_km.csv and load_bedmap_catalog fills line_km from it instead of catalog geometry. - plot_survey_density.load_bedmap uses the shared functions. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_019UvrC4evHHkT7qWWL96sD5 --- README.md | 11 +++- bedmap_campaign_km.csv | 114 ++++++++++++++++++++++++++++++++++ bedmap_common.py | 85 ++++++++++++++++++++++++- plot_bedmap_availability.py | 5 +- plot_combined_availability.py | 1 - plot_survey_density.py | 31 +++------ 6 files changed, 214 insertions(+), 33 deletions(-) create mode 100644 bedmap_campaign_km.csv diff --git a/README.md b/README.md index 37b8348..9365e0c 100644 --- a/README.md +++ b/README.md @@ -20,9 +20,14 @@ Or edit `pyproject.toml` to point `xopr` at your local path. All BedMap consumers share `bedmap_common.py`, which applies the campaign exclusion list, drops BM2 files superseded by renamed BM3 files, keeps one -version of campaigns present in both BedMap2 and BedMap3, and parses temporal -metadata (a missing or year-9999 end date falls back to the start date). -Change those rules there, not in individual plots. +version of campaigns present in both BedMap2 and BedMap3, parses temporal +metadata (a missing or year-9999 end date falls back to the start date), and +defines the gap rule (`MAX_POINT_SPACING_M`, 1.5 km) that decides which +consecutive points form a flight line. Per-campaign line-km for the bar +charts come from that rule applied to the point data, stored in +`bedmap_campaign_km.csv` (rebuild with `uv run python bedmap_common.py` +after changing the rules; needs the point cache), not from the simplified +catalog geometry. Change those rules there, not in individual plots. Direct provider sources live in `extra_sources.py` and are added to every BedMap-based figure (bar charts, survey density, coastal gaps; `--no-extra` diff --git a/bedmap_campaign_km.csv b/bedmap_campaign_km.csv new file mode 100644 index 0000000..bb8108e --- /dev/null +++ b/bedmap_campaign_km.csv @@ -0,0 +1,114 @@ +name,line_km +BAS_1994_Evans_AIR_BM2,12179.262650741326 +BAS_1998_Dufek_AIR_BM2,4149.363483508933 +BAS_2001_Bailey-Slessor_AIR_BM2,9417.102981384793 +BAS_2001_MAMOG_AIR_BM2,2120.9609719280384 +BAS_2001_TORUS_AIR_BM2,7172.650540529199 +BAS_2004_BBAS_AIR_BM2,34282.60491095511 +BAS_2005_WISE-ISODYN_AIR_BM2,59688.18212625556 +BAS_2006_GRADES-IMAGE_AIR_BM2,22864.09683745935 +BAS_2007_AGAP_AIR_BM2,53245.22499105713 +BAS_2007_Lake-Ellsworth_GRN_BM3,291.39961512183146 +BAS_2007_Rutford_GRN_BM3,1274.9219945992415 +BAS_2007_TIGRIS_GRN_BM2,1105.800692939278 +BAS_2008_Lake-Ellsworth_GRN_BM3,545.0559380028765 +BAS_2009_FERRIGNO_GRN_BM2,1246.2029703357766 +BAS_2010_IMAFI_AIR_BM3,21935.345429941688 +BAS_2010_PIG_AIR_BM2,1420.557094373895 +BAS_2011_Adelaide_AIR_BM3,2496.9405813433536 +BAS_2012_Castle_GRN_BM3,212.03544768235724 +BAS_2012_ICEGRAV_AIR_BM3,25010.898751509045 +BAS_2013_ISTAR_GRN_BM3,2558.8971836098713 +BAS_2015_FISS_AIR_BM3,5496.2353004141905 +BAS_2015_POLARGAP_AIR_BM3,28924.060046695704 +BAS_2016_FISS_AIR_BM3,14099.682974756211 +BAS_2017_English-Coast_AIR_BM3,3719.7530757162654 +BAS_2018_Thwaites_AIR_BM3,4254.490633743226 +BAS_2019_Thwaites_AIR_BM3,4496.8047487926015 +BEDMAP1_1966-2000_AIR_BM1,512751.4748787665 +BGR_1999_GANOVEX-VIII-Matusevich_AIR_BM2,6783.0495726106765 +BGR_1999_GANOVEX-VIII-Mertz_AIR_BM2,1300.4178330820591 +BGR_2002_PCMEGA_AIR_BM2,31769.13425950711 +CECS_2006_Subglacial-Lake-CECs_GRN_BM3,1939.2119516739808 +CRESIS_2013_Siple-Coast_AIR_BM3,17274.972127928067 +INGV_1997_ITASE_AIR_BM2,3256.8442284470516 +INGV_1997_Talos-Dome_AIR_BM3,2263.008632785903 +INGV_1999_Talos-Dome_AIR_BM3,11153.019749818553 +INGV_2001_Talos-Dome_AIR_BM3,3817.8735007922646 +INGV_2003_Talos-Dome_AIR_BM3,22061.592422766287 +KOPRI_2017_KRT1_AIR_BM3,4339.559925100045 +KOPRI_2018_KRT2_AIR_BM3,1517.4317778335753 +LDEO_2007_AGAP-GAMBIT_AIR_BM2,32609.428173537333 +LDEO_2007_Recovery-Lakes_AIR_BM2,2877.7324600031016 +LDEO_2015_ROSETTA_AIR_BM3,43523.444895789056 +NASA_2002_ICEBRIDGE_AIR_BM2,6238.885501372943 +NASA_2004_ICEBRIDGE_AIR_BM2,4290.652034216112 +NASA_2009_ICEBRIDGE_AIR_BM2,25425.4062116988 +NASA_2010_ICEBRIDGE_AIR_BM2,20071.19778493137 +NASA_2011_ICEBRIDGE_AIR_BM2,65228.757506913884 +NASA_2012_ICEBRIDGE_AIR_BM2,32623.679849093303 +NASA_2013_ICEBRIDGE_AIR_BM3,14166.211692343495 +NASA_2014_ICEBRIDGE_AIR_BM3,52463.17606975831 +NASA_2016_ICEBRIDGE_AIR_BM3,55826.74931727524 +NASA_2017_ICEBRIDGE_AIR_BM3,30979.401647501687 +NASA_2018_ICEBRIDGE_AIR_BM3,40577.38595892682 +NASA_2019_ICEBRIDGE_AIR_BM3,20930.313516593418 +NIPR_1992_JARE33_GRN_BM3,373.2783262633138 +NIPR_1996_JARE37_GRN_BM3,1878.6244437050004 +NIPR_1999_JARE40_GRN_BM3,172.83667006028256 +NIPR_2007_JARE49_GRN_BM3,1438.3714194173222 +NIPR_2007_JASE_GRN_BM3,2300.734965231171 +NIPR_2012_JARE54_GRN_BM3,126.65998551904391 +NIPR_2017_JARE59_GRN_BM3,2014.7194283145916 +NIPR_2018_JARE60_GRN_BM3,8905.314937862864 +NPI_2008_BELISSIMA_GRN_BM2,66.30617524699694 +NPI_2010_SRM_AIR_BM2,3270.4030055834273 +NPI_2012_ICERISES_GRN_BM3,327.41267143864115 +NPI_2015_POLARGAP_AIR_BM3,2815.3667153382175 +NPI_2016_MADICE_GRN_BM3,610.904432433296 +PRIC_2004_CHINARE-21_GRN_BM2,2559.982587907926 +PRIC_2007_CHINARE-24_GRN_BM2,603.3034937502166 +PRIC_2015_CHA1_AIR_BM3,40177.03816930688 +PRIC_2016_CHA2_AIR_BM3,21004.283508419878 +PRIC_2017_CHA3_AIR_BM3,25446.214595007812 +PRIC_2018_CHA4_AIR_BM3,15052.899018330012 +RNRF_1975_Filchner-Ronne_SEI_BM3,1.2835160741144491 +RNRF_1975_Lazarev_SEI_BM3,9.422076160352173 +RNRF_2003_AMSap5_AIR_BM3,4955.185928409161 +RNRF_2004_AMSap5_AIR_BM3,6360.120969857368 +RNRF_2004_Mirny-Vostok_AIR_BM3,1415.5937777432682 +RNRF_2005_AMSap5_AIR_BM3,4247.575813578208 +RNRF_2006_Komsom-Vostok_AIR_BM3,481.9184365624733 +RNRF_2006_RAEap5_AIR_BM3,6001.475438810396 +RNRF_2007_AMSap5_AIR_BM3,5729.24149259773 +RNRF_2008_AMSap5_AIR_BM3,5151.903848739583 +RNRF_2009_RAEap5_AIR_BM3,5002.741833256816 +RNRF_2010_RAE_AIR_BM3,7363.233995436161 +RNRF_2011_RAE_AIR_BM3,6977.69250118063 +RNRF_2013_RAE_AIR_BM3,14361.90076059565 +RNRF_2013_Vostok-Progress_AIR_BM2,1226.6973135938579 +RNRF_2014_RAE_AIR_BM3,6450.733145429845 +RNRF_2015_RAE_AIR_BM3,6096.79731173056 +RNRF_2016_RAE_AIR_BM3,6490.6332234855 +RNRF_2017_RAE_AIR_BM3,6839.40239536566 +RNRF_2018_RAE_AIR_BM3,5864.027618573344 +RNRF_2019_RAE_AIR_BM3,5654.696345256103 +STANFORD_1971_SPRI-NSF-TUD_AIR_BM3,7286.122531415287 +STOLAF_1994_Siple-Dome_GRN_BM2,143.5455410337653 +STOLAF_2001_ITASE-Byrd-Ellsworth_GRN_BM2,879.0093416878482 +STOLAF_2001_ITASE-Ellsworth_GRN_BM2,145.10345224251216 +STOLAF_2002_ITASE-Byrd-South-Pole_GRN_BM2,897.8164427905398 +STOLAF_2002_ITASE-Hercules-Dome_GRN_BM2,121.7829483582652 +UCANTERBURY_2008_Darwin-Hatherton_GRN_BM2,216.05856338281825 +ULB_2012_BEWISE_GRN_BM3,38.48529497266142 +ULB_2012_ICECON_GRN_BM3,288.4447429421593 +UTIG_1991_CASERTZ_AIR_BM2,29808.045520122367 +UTIG_1998_West-Marie-Byrd-Land_AIR_BM2,51534.00051425103 +UTIG_1999_SOAR-LVS-WLK_AIR_BM2,152837.89738200157 +UTIG_2000_Robb-Glacier_AIR_BM2,832.2557513557491 +UTIG_2009_Darwin-Hatherton_AIR_BM3,698.5972199718469 +UTIG_2010_ICECAP_AIR_BM3,201552.18361124926 +UTIG_2013_GIMBLE_AIR_BM3,11562.400626045273 +UTIG_2015_EAGLE_AIR_BM3,21958.623083374136 +UTIG_2016_OLDICE_AIR_BM3,9142.624259238333 +UWASHINGTON_2018_South-Pole-Lake_GRN_BM3,187.95916516155825 diff --git a/bedmap_common.py b/bedmap_common.py index 1f142df..723abf1 100644 --- a/bedmap_common.py +++ b/bedmap_common.py @@ -6,12 +6,22 @@ """ import re +from pathlib import Path import numpy as np import pandas as pd -from pyproj import Geod +from pyproj import Geod, Transformer from xopr.bedmap.query import query_bedmap_catalog +# Consecutive BedMap points further apart than this are not joined into a +# flight line. Applied to the density grids and to the per-campaign line-km +# used by the bar charts. 1.5 km keeps the airborne surveys sampled near 1 km +# (SOAR-LVS-WLK, CASERTZ, MAMOG: 95th-percentile spacing <= 1.35 km) and still +# drops ground traverses and the 25 km-spaced GANOVEX files (p95 >= 6.5 km); +# see claude_notes/bedmap_gap_rule_options.csv. +MAX_POINT_SPACING_M = 1500 +CAMPAIGN_KM_FILE = Path(__file__).parent / "bedmap_campaign_km.csv" + # Campaigns dropped everywhere. CRESIS_2009_Thwaites is a MUSIC swath product # (~30 cross-track bed picks per radar record) and the Vostok file is scattered # points; neither can be chained into flight lines. @@ -60,8 +70,11 @@ def load_bedmap_catalog(collections=("bedmap2", "bedmap3"), exclude=True): metadata into ``ts``/``te``. A missing or sentinel end date (e.g. year 9999) falls back to the start date. + ``line_km`` is the gap-filtered length of the campaign's point data + (see ``campaign_point_km``), not the simplified catalog geometry. + Columns: collection, geometry, name, base_name, institution, - temporal_start, temporal_end, ts, te. + temporal_start, temporal_end, ts, te, line_km. """ cat = query_bedmap_catalog(collections=list(collections)) props = pd.DataFrame(cat["properties"].tolist(), index=cat.index) @@ -81,9 +94,60 @@ def load_bedmap_catalog(collections=("bedmap2", "bedmap3"), exclude=True): te = pd.to_datetime(df["temporal_end"], format="ISO8601", errors="coerce") bad_end = te.isna() | (te.dt.year > 2100) df["te"] = te.where(~bad_end, df["ts"]) + km = campaign_point_km() + missing = df.loc[~df["name"].isin(km.index), "name"].tolist() + if missing: + print(f"bedmap_common: no point-based km for {missing}; using catalog geometry") + df["line_km"] = df["name"].map(km).fillna(df["geometry"].apply(geod_km)) return df.reset_index(drop=True) +def load_bedmap_points(collections=("bedmap1", "bedmap2", "bedmap3"), local_cache=True): + """Kept BedMap points (lon, lat, source_file, row) sorted along track.""" + from xopr.bedmap import fetch_bedmap, query_bedmap + fetch_bedmap() + df = query_bedmap(collections=list(collections), + columns=["lon", "lat", "source_file", "row"], + local_cache=local_cache, show_progress=True) + n_files = df["source_file"].nunique() + df = df[df["source_file"].isin(kept_source_files(collections))] + print(f" {len(df)} BedMap points from {df['source_file'].nunique()} files " + f"({n_files - df['source_file'].nunique()} excluded or duplicate files dropped)") + return df.sort_values(["source_file", "row"]) + + +def point_segments(df, epsg="EPSG:3031", max_point_spacing_m=MAX_POINT_SPACING_M): + """Segments between consecutive points of one file, shorter than the gap limit. + + Returns x1, y1, x2, y2 in ``epsg`` plus the segment lengths in geodesic + km and the source_file of each segment. + """ + tf = Transformer.from_crs("EPSG:4326", epsg, always_xy=True) + lon, lat = df["lon"].to_numpy(), df["lat"].to_numpy() + xs, ys = tf.transform(lon, lat) + same_file = df["source_file"].to_numpy()[:-1] == df["source_file"].to_numpy()[1:] + dists = np.hypot(np.diff(xs), np.diff(ys)) + ok = same_file & (dists < max_point_spacing_m) & (dists > 0) + km = _GEOD.inv(lon[:-1][ok], lat[:-1][ok], lon[1:][ok], lat[1:][ok])[2] / 1000.0 + return xs[:-1][ok], ys[:-1][ok], xs[1:][ok], ys[1:][ok], km, df["source_file"].to_numpy()[:-1][ok] + + +def campaign_point_km(rebuild=False): + """Gap-filtered line-km per BedMap file, from ``bedmap_campaign_km.csv``. + + The table is rebuilt from the point data (all three BedMap versions, + ~90M points, needs the point cache) when missing or ``rebuild`` is set; + run ``python bedmap_common.py`` to refresh it after changing the rules. + """ + if CAMPAIGN_KM_FILE.exists() and not rebuild: + return pd.read_csv(CAMPAIGN_KM_FILE, index_col="name")["line_km"] + df = load_bedmap_points() + *_, seg_km, files = point_segments(df) + km = pd.Series(seg_km).groupby(files).sum().rename("line_km").rename_axis("name") + km.to_csv(CAMPAIGN_KM_FILE) + return km + + def dropped_names(): """BedMap file names removed everywhere (excluded, superseded, substituted).""" subst = {n for names in OPR_SUBSTITUTIONS.values() for n in names} @@ -108,7 +172,17 @@ def kept_source_files(collections=("bedmap1", "bedmap2", "bedmap3")): ``source_file`` in the point parquet files matches the catalog ``name``. """ - return set(load_bedmap_catalog(collections)["name"]) + return set(_catalog_names(collections)) + + +def _catalog_names(collections): + """Kept catalog names (exclusions + version dedup) without the km table.""" + cat = query_bedmap_catalog(collections=list(collections)) + df = pd.DataFrame({"name": [p["name"] for p in cat["properties"]]}) + df = df[~df["name"].isin(dropped_names()) + & ~df["name"].str.startswith(EXTERNAL_SOURCE_PREFIXES)] + df["base"] = df["name"].str.replace(_VERSION_SUFFIX, "", regex=True) + return df.sort_values("name").drop_duplicates("base", keep="last")["name"] def geod_km(geometry): @@ -121,3 +195,8 @@ def geod_km(geometry): if len(c) > 1: total += _GEOD.inv(c[:-1, 0], c[:-1, 1], c[1:, 0], c[1:, 1])[2].sum() return total / 1000.0 + + +if __name__ == "__main__": + km = campaign_point_km(rebuild=True) + print(f"wrote {CAMPAIGN_KM_FILE}: {len(km)} files, {km.sum():,.0f} km") diff --git a/plot_bedmap_availability.py b/plot_bedmap_availability.py index 1592b4e..c6d36a7 100644 --- a/plot_bedmap_availability.py +++ b/plot_bedmap_availability.py @@ -12,7 +12,7 @@ import matplotlib.pyplot as plt import pandas as pd -from bedmap_common import campaign_years, geod_km, load_bedmap_catalog +from bedmap_common import campaign_years, load_bedmap_catalog from extra_sources import load_extra_campaigns SCRIPT_DIR = Path(__file__).parent @@ -42,9 +42,8 @@ def institution_to_country(name): # Query bedmap2 and bedmap3 catalogs (skip bedmap1, matching reference). -# Exclusions, BM2/BM3 dedup and date parsing live in bedmap_common. +# Exclusions, dedup, date parsing and gap-filtered line-km live in bedmap_common. df = load_bedmap_catalog(["bedmap2", "bedmap3"]) -df["line_km"] = df["geometry"].apply(geod_km) df = pd.concat([df, load_extra_campaigns()], ignore_index=True) # Distribute every campaign evenly across its years (see campaign_years) diff --git a/plot_combined_availability.py b/plot_combined_availability.py index 369a1b3..192933d 100644 --- a/plot_combined_availability.py +++ b/plot_combined_availability.py @@ -38,7 +38,6 @@ print("Querying BedMap catalogs...") # Exclusions, BM2/BM3 dedup and date parsing live in bedmap_common. bm = load_bedmap_catalog(["bedmap2", "bedmap3"]) - bm["line_km"] = bm["geometry"].apply(geod_km) bm = pd.concat([bm, load_extra_campaigns()], ignore_index=True) bm_rows = [] diff --git a/plot_survey_density.py b/plot_survey_density.py index fb08ebc..9b997a0 100644 --- a/plot_survey_density.py +++ b/plot_survey_density.py @@ -67,15 +67,16 @@ } -def load_bedmap(epsg, max_point_spacing_m=1000, local_cache=True, include_extra=True): - """Load BedMap point data (plus extra_sources tracks) as projected segment endpoints. +def load_bedmap(epsg, local_cache=True, include_extra=True): + """BedMap point data (plus extra_sources tracks) as projected segment endpoints. + + Uses the shared gap rule in ``bedmap_common`` (``MAX_POINT_SPACING_M``), so + the grids and the bar charts count the same line-km. Parameters ---------- epsg : str Target CRS (e.g. "EPSG:3031"). - max_point_spacing_m : float - Discard segments longer than this (gap between survey points). include_extra : bool Append the direct-source tracks from extra_sources (AWI, KOPRI, xOPR substitutes), which replace the corresponding BedMap files. @@ -85,25 +86,9 @@ def load_bedmap(epsg, max_point_spacing_m=1000, local_cache=True, include_extra= x1, y1, x2, y2 : np.ndarray Segment endpoint arrays in the target CRS. """ - from xopr.bedmap import query_bedmap, fetch_bedmap - tf = Transformer.from_crs("EPSG:4326", epsg, always_xy=True) - fetch_bedmap() - df = query_bedmap(collections=["bedmap1", "bedmap2", "bedmap3"], - columns=["lon", "lat", "source_file", "row"], - local_cache=local_cache, show_progress=True) - # Drop excluded campaigns and BM2/BM3 duplicates (see bedmap_common) - from bedmap_common import kept_source_files - n_files = df["source_file"].nunique() - df = df[df["source_file"].isin(kept_source_files(["bedmap1", "bedmap2", "bedmap3"]))] - print(f" {len(df)} BedMap points from {df['source_file'].nunique()} files " - f"({n_files - df['source_file'].nunique()} excluded or duplicate files dropped)") - df = df.sort_values(["source_file", "row"]) - xs, ys = tf.transform(df["lon"].values, df["lat"].values) - # Mask transitions between files so we don't connect unrelated points - same_file = df["source_file"].values[:-1] == df["source_file"].values[1:] - dists = np.sqrt(np.diff(xs)**2 + np.diff(ys)**2) - ok = same_file & (dists < max_point_spacing_m) & (dists > 0) - segs = [xs[:-1][ok], ys[:-1][ok], xs[1:][ok], ys[1:][ok]] + from bedmap_common import load_bedmap_points, point_segments + df = load_bedmap_points(local_cache=local_cache) + segs = list(point_segments(df, epsg)[:4]) if include_extra: from extra_sources import extra_segments ex = extra_segments(epsg) From a79ab55c4c3fed763b8339e2412228df7f85f655 Mon Sep 17 00:00:00 2001 From: Thomas Teisberg Date: Tue, 15 Sep 2026 15:34:37 -0700 Subject: [PATCH 04/10] Country chart: --no-haps variant, y axis in thousands of km, larger fonts Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_019UvrC4evHHkT7qWWL96sD5 --- README.md | 2 ++ plot_bedmap_availability.py | 41 ++++++++++++++++++++++--------------- 2 files changed, 26 insertions(+), 17 deletions(-) diff --git a/README.md b/README.md index 9365e0c..098a9b5 100644 --- a/README.md +++ b/README.md @@ -90,6 +90,8 @@ Stacked bar chart of BedMap line-km per year, colored by country of origin. ```bash uv run python plot_bedmap_availability.py # -> outputs/bedmap_data_availability.png +uv run python plot_bedmap_availability.py --no-haps +# -> outputs/bedmap_data_availability_nohaps.png (without the HAPS capability line) ``` ### `plot_opr_availability.py` diff --git a/plot_bedmap_availability.py b/plot_bedmap_availability.py index c6d36a7..1c6bc0e 100644 --- a/plot_bedmap_availability.py +++ b/plot_bedmap_availability.py @@ -7,6 +7,7 @@ temporal_start..temporal_end calendar years; dated sources go to their season. """ +import argparse from pathlib import Path import matplotlib.pyplot as plt @@ -15,6 +16,10 @@ from bedmap_common import campaign_years, load_bedmap_catalog from extra_sources import load_extra_campaigns +p = argparse.ArgumentParser(description=__doc__) +p.add_argument("--no-haps", action="store_true", help="Omit the HAPS capability line") +args = p.parse_args() + SCRIPT_DIR = Path(__file__).parent OUT_DIR = SCRIPT_DIR / "outputs" OUT_DIR.mkdir(exist_ok=True) @@ -64,26 +69,28 @@ def institution_to_country(name): pivot = pivot.reindex(range(2000, 2026), fill_value=0) pivot = pivot.reindex(columns=[c for c in COUNTRY_ORDER if c in pivot.columns]) -# Plot +# Plot (thousands of km) fig, ax = plt.subplots(figsize=(14, 7)) -pivot.plot.bar(stacked=True, ax=ax, - color=[COUNTRY_COLORS[c] for c in pivot.columns], - width=0.8, edgecolor="none") - -# HAPS reference line -haps_km = 153000 -ax.axhline(haps_km, color="red", linestyle="--", linewidth=2) -ax.annotate("Capability of 1 HAPS UAV, 11 week mission", - xy=(0.35, haps_km + 2000), xycoords=("axes fraction", "data"), - fontsize=14, color="red", fontweight="bold", ha="center") - -ax.set_title("Line-km of global Antarctic airborne radar surveying", fontsize=16) -ax.set_xlabel("year", fontsize=13) -ax.set_ylabel("IPR surveying flight kilometers", fontsize=13) -ax.legend(title="Country", fontsize=11, title_fontsize=12, +(pivot / 1000).plot.bar(stacked=True, ax=ax, + color=[COUNTRY_COLORS[c] for c in pivot.columns], + width=0.8, edgecolor="none") + +if not args.no_haps: + haps_km = 153 + ax.axhline(haps_km, color="red", linestyle="--", linewidth=2) + ax.annotate("Capability of 1 HAPS UAV, 11 week mission", + xy=(0.35, haps_km + 2), xycoords=("axes fraction", "data"), + fontsize=16, color="red", fontweight="bold", ha="center") + +ax.set_title("Line-km of global Antarctic airborne radar surveying", fontsize=20) +ax.set_xlabel("year", fontsize=17) +ax.set_ylabel("IPR surveying flight kilometers (thousands)", fontsize=17) +ax.tick_params(axis="both", labelsize=14) +ax.legend(title="Country", fontsize=14, title_fontsize=15, loc="upper right", framealpha=0.9) plt.tight_layout() -out_path = OUT_DIR / "bedmap_data_availability.png" +suffix = "_nohaps" if args.no_haps else "" +out_path = OUT_DIR / f"bedmap_data_availability{suffix}.png" plt.savefig(out_path, dpi=150, bbox_inches="tight") plt.show() print(f"Saved to {out_path}") From 4c9957b8d54586b4602e65e97aafeac098b31a23 Mon Sep 17 00:00:00 2001 From: Thomas Teisberg Date: Tue, 15 Sep 2026 15:51:08 -0700 Subject: [PATCH 05/10] Country chart: stack USA at the bottom Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_019UvrC4evHHkT7qWWL96sD5 --- plot_bedmap_availability.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/plot_bedmap_availability.py b/plot_bedmap_availability.py index 1c6bc0e..7aa584d 100644 --- a/plot_bedmap_availability.py +++ b/plot_bedmap_availability.py @@ -37,7 +37,7 @@ "USA": "tab:blue", "UK": "tab:green", "Germany": "tab:orange", "Russia": "tab:red", "China": "tab:purple", "Other": "gray", } -COUNTRY_ORDER = ["Other", "China", "Russia", "Germany", "UK", "USA"] +COUNTRY_ORDER = ["USA", "UK", "Germany", "Russia", "China", "Other"] # bottom to top def institution_to_country(name): From 1aa4510a64dd6e22ae9f8b8aa1efc52d3fafb99b Mon Sep 17 00:00:00 2001 From: Thomas Teisberg Date: Tue, 15 Sep 2026 16:05:30 -0700 Subject: [PATCH 06/10] Spread the SOAR 1996-2001 BedMap file over its real years (TEMPORAL_OVERRIDES) Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_019UvrC4evHHkT7qWWL96sD5 --- README.md | 3 ++- bedmap_common.py | 12 +++++++++++- 2 files changed, 13 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 098a9b5..0d8aec3 100644 --- a/README.md +++ b/README.md @@ -21,7 +21,8 @@ Or edit `pyproject.toml` to point `xopr` at your local path. All BedMap consumers share `bedmap_common.py`, which applies the campaign exclusion list, drops BM2 files superseded by renamed BM3 files, keeps one version of campaigns present in both BedMap2 and BedMap3, parses temporal -metadata (a missing or year-9999 end date falls back to the start date), and +metadata (a missing or year-9999 end date falls back to the start date; +`TEMPORAL_OVERRIDES` fixes known bad ranges such as the SOAR 1996-2001 file), and defines the gap rule (`MAX_POINT_SPACING_M`, 1.5 km) that decides which consecutive points form a flight line. Per-campaign line-km for the bar charts come from that rule applied to the point data, stored in diff --git a/bedmap_common.py b/bedmap_common.py index 723abf1..66bbb0e 100644 --- a/bedmap_common.py +++ b/bedmap_common.py @@ -57,6 +57,13 @@ # instead of their BedMap submissions. EXTERNAL_SOURCE_PREFIXES = ("AWI_",) +# Files whose catalog dates are known to be wrong: name -> (start, end). +# UTIG_1999_SOAR-LVS-WLK bundles the SOAR surveys of 1996-2001 (TAM/PPT, +# WLK, BSB, IRE, WAZ/TKD, LVS; USAP-DC 601588) under a 1999-2000 placeholder. +TEMPORAL_OVERRIDES = { + "UTIG_1999_SOAR-LVS-WLK_AIR_BM2": ("1996-01-01", "2001-12-31"), +} + _VERSION_SUFFIX = re.compile(r"_BM[123]$") _GEOD = Geod(ellps="WGS84") @@ -68,7 +75,8 @@ def load_bedmap_catalog(collections=("bedmap2", "bedmap3"), exclude=True): ``EXTERNAL_SOURCE_PREFIXES``, keeps the newest BedMap version when a campaign appears in several (BM3 over BM2 over BM1), and parses temporal metadata into ``ts``/``te``. A missing or sentinel end date (e.g. year - 9999) falls back to the start date. + 9999) falls back to the start date; ``TEMPORAL_OVERRIDES`` replaces known + bad ranges. ``line_km`` is the gap-filtered length of the campaign's point data (see ``campaign_point_km``), not the simplified catalog geometry. @@ -94,6 +102,8 @@ def load_bedmap_catalog(collections=("bedmap2", "bedmap3"), exclude=True): te = pd.to_datetime(df["temporal_end"], format="ISO8601", errors="coerce") bad_end = te.isna() | (te.dt.year > 2100) df["te"] = te.where(~bad_end, df["ts"]) + for name, (start, end) in TEMPORAL_OVERRIDES.items(): + df.loc[df["name"] == name, ["ts", "te"]] = [pd.Timestamp(start, tz="UTC"), pd.Timestamp(end, tz="UTC")] km = campaign_point_km() missing = df.loc[~df["name"].isin(km.index), "name"].tolist() if missing: From ff5d4494f497eb7291dd380c8f100c342b8e8abe Mon Sep 17 00:00:00 2001 From: Thomas Teisberg Date: Wed, 16 Sep 2026 11:59:02 -0700 Subject: [PATCH 07/10] Pages: add the by-country chart (no HAPS line) Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_019UvrC4evHHkT7qWWL96sD5 --- .github/workflows/generate-figures.yml | 3 +++ index.html | 10 +++++++++- 2 files changed, 12 insertions(+), 1 deletion(-) diff --git a/.github/workflows/generate-figures.yml b/.github/workflows/generate-figures.yml index b686e93..cb0441a 100644 --- a/.github/workflows/generate-figures.yml +++ b/.github/workflows/generate-figures.yml @@ -32,6 +32,9 @@ jobs: - name: Generate combined availability figure run: uv run python plot_combined_availability.py + - name: Generate by-country availability figure + run: uv run python plot_bedmap_availability.py --no-haps + - name: Generate BedMap Antarctica survey density run: uv run python plot_survey_density.py --source bedmap --region antarctica diff --git a/index.html b/index.html index 8ee8b61..868054b 100644 --- a/index.html +++ b/index.html @@ -20,7 +20,8 @@

Radar Data Availability Figures

These are some useful figures showing airborne radar sounder data availability and survey density for the Antarctic and Greenland ice sheets, using data from - xOPR and BedMap catalogs. + xOPR and BedMap catalogs, supplemented with + flight tracks published directly by AWI and the KOPRI/UTIG helicopter surveys. Generated automatically on each push to main.

@@ -35,6 +36,13 @@

Combined Data Availability

uv run python plot_combined_availability.py
+

Antarctic Airborne Radar Surveying by Country

+
+ Antarctic airborne radar line-km per year by country +
Line-km of Antarctic airborne radar sounding per year, stacked by the country of the collecting institution. BedMap campaigns are spread over their catalog date range; the AWI, KOPRI and xOPR-sourced tracks are placed in their field season.
+
+ 
uv run python plot_bedmap_availability.py --no-haps
+

Survey Density: BedMap Antarctica

BedMap Antarctica survey density From bc21cda3beab879264275ee76697ef837688fb88 Mon Sep 17 00:00:00 2001 From: Thomas Teisberg Date: Wed, 16 Sep 2026 14:08:06 -0700 Subject: [PATCH 08/10] Country chart: opt-in --reference-line KM / --reference-label instead of --no-haps Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_019UvrC4evHHkT7qWWL96sD5 --- .github/workflows/generate-figures.yml | 2 +- README.md | 4 ++-- index.html | 4 ++-- plot_bedmap_availability.py | 21 +++++++++++++-------- 4 files changed, 18 insertions(+), 13 deletions(-) diff --git a/.github/workflows/generate-figures.yml b/.github/workflows/generate-figures.yml index cb0441a..1736035 100644 --- a/.github/workflows/generate-figures.yml +++ b/.github/workflows/generate-figures.yml @@ -33,7 +33,7 @@ jobs: run: uv run python plot_combined_availability.py - name: Generate by-country availability figure - run: uv run python plot_bedmap_availability.py --no-haps + run: uv run python plot_bedmap_availability.py - name: Generate BedMap Antarctica survey density run: uv run python plot_survey_density.py --source bedmap --region antarctica diff --git a/README.md b/README.md index 0d8aec3..0e667ef 100644 --- a/README.md +++ b/README.md @@ -91,8 +91,8 @@ Stacked bar chart of BedMap line-km per year, colored by country of origin. ```bash uv run python plot_bedmap_availability.py # -> outputs/bedmap_data_availability.png -uv run python plot_bedmap_availability.py --no-haps -# -> outputs/bedmap_data_availability_nohaps.png (without the HAPS capability line) +uv run python plot_bedmap_availability.py --reference-line 153000 --reference-label "Capability of 1 HAPS UAV, 11 week mission" +# -> outputs/bedmap_data_availability_ref.png (with a labelled reference line) ``` ### `plot_opr_availability.py` diff --git a/index.html b/index.html index 868054b..eb8c8a9 100644 --- a/index.html +++ b/index.html @@ -38,10 +38,10 @@

Combined Data Availability

Antarctic Airborne Radar Surveying by Country

- Antarctic airborne radar line-km per year by country + Antarctic airborne radar line-km per year by country
Line-km of Antarctic airborne radar sounding per year, stacked by the country of the collecting institution. BedMap campaigns are spread over their catalog date range; the AWI, KOPRI and xOPR-sourced tracks are placed in their field season.
- 
uv run python plot_bedmap_availability.py --no-haps
+ 
uv run python plot_bedmap_availability.py

Survey Density: BedMap Antarctica

diff --git a/plot_bedmap_availability.py b/plot_bedmap_availability.py index 7aa584d..845b525 100644 --- a/plot_bedmap_availability.py +++ b/plot_bedmap_availability.py @@ -17,7 +17,11 @@ from extra_sources import load_extra_campaigns p = argparse.ArgumentParser(description=__doc__) -p.add_argument("--no-haps", action="store_true", help="Omit the HAPS capability line") +p.add_argument("--reference-line", type=float, metavar="KM", + help="Draw a horizontal reference line at this many line-km " + "(e.g. 153000 for one HAPS UAV over an 11 week mission)") +p.add_argument("--reference-label", default="", + help="Label for the reference line") args = p.parse_args() SCRIPT_DIR = Path(__file__).parent @@ -75,12 +79,13 @@ def institution_to_country(name): color=[COUNTRY_COLORS[c] for c in pivot.columns], width=0.8, edgecolor="none") -if not args.no_haps: - haps_km = 153 - ax.axhline(haps_km, color="red", linestyle="--", linewidth=2) - ax.annotate("Capability of 1 HAPS UAV, 11 week mission", - xy=(0.35, haps_km + 2), xycoords=("axes fraction", "data"), - fontsize=16, color="red", fontweight="bold", ha="center") +if args.reference_line is not None: + ref = args.reference_line / 1000 + ax.axhline(ref, color="red", linestyle="--", linewidth=2) + if args.reference_label: + ax.annotate(args.reference_label, xy=(0.35, ref + 2), + xycoords=("axes fraction", "data"), fontsize=16, + color="red", fontweight="bold", ha="center") ax.set_title("Line-km of global Antarctic airborne radar surveying", fontsize=20) ax.set_xlabel("year", fontsize=17) @@ -89,7 +94,7 @@ def institution_to_country(name): ax.legend(title="Country", fontsize=14, title_fontsize=15, loc="upper right", framealpha=0.9) plt.tight_layout() -suffix = "_nohaps" if args.no_haps else "" +suffix = "_ref" if args.reference_line is not None else "" out_path = OUT_DIR / f"bedmap_data_availability{suffix}.png" plt.savefig(out_path, dpi=150, bbox_inches="tight") plt.show() From 6e4a8f9eba76463887a0d424b3696dcc88790972 Mon Sep 17 00:00:00 2001 From: Thomas Teisberg Date: Wed, 16 Sep 2026 14:24:15 -0700 Subject: [PATCH 09/10] Fetch AWI tracks in filtered pages with retries (CI download was truncated) The single 245 MB GetFeature response was cut off on the CI runner. Use a server-side radar_system filter (8.8k of 40.7k features) and 2000-feature pages, and retry truncated downloads for all sources. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_019UvrC4evHHkT7qWWL96sD5 --- extra_sources.py | 44 +++++++++++++++++++++++++++++++++++++------- opr_tracks.py | 7 ++----- 2 files changed, 39 insertions(+), 12 deletions(-) diff --git a/extra_sources.py b/extra_sources.py index ee89a60..4132654 100644 --- a/extra_sources.py +++ b/extra_sources.py @@ -13,6 +13,8 @@ Raw downloads are cached under ``radar_cache/`` (gitignored). """ +import time +import urllib.parse import urllib.request from pathlib import Path @@ -24,9 +26,9 @@ from bedmap_common import geod_km CACHE = Path(__file__).parent / "radar_cache" -AWI_WFS = ("https://maps.awi.de/services/common/radar/ows?service=WFS&version=2.0.0" - "&request=GetFeature&typeNames=radar:tracks_3031&outputFormat=application/json") +AWI_WFS = "https://maps.awi.de/services/common/radar/ows" AWI_SOUNDERS = {"EMR", "UWB", "UWBM"} # ACCU/SNOW/ASIRAS are accumulation radars +AWI_PAGE = 2000 # features per WFS request; the full layer in one response gets cut off AWI_GENERIC_DOI = "10.1594/PANGAEA.972094" # collection bibliography, not a dataset release QICERADAR_GPKG = "https://zenodo.org/api/records/21964546/files/qiceradar_antarctic_index.gpkg/content" # KOPRI Araon cruises are annual and visit Thwaites every other year (last @@ -34,10 +36,21 @@ KOPRI_SEASONS = {"ASE2": 2017, "ASE3": 2019, "ASE4": 2021, "ASE5": 2023, "ASE6": 2025} -def _fetch(url, path): +def _fetch(url, path, attempts=4): + """Download ``url`` to ``path`` unless present; retry on truncated responses.""" if not path.exists(): path.parent.mkdir(parents=True, exist_ok=True) - urllib.request.urlretrieve(url, path) + tmp = path.with_suffix(path.suffix + ".part") + for i in range(attempts): + try: + urllib.request.urlretrieve(url, tmp) + tmp.rename(path) + break + except Exception as e: # noqa: BLE001 - network errors of any kind + if i == attempts - 1: + raise + print(f" retrying {url.split('?')[0]} after {type(e).__name__}") + time.sleep(5 * (i + 1)) return path @@ -46,11 +59,28 @@ def _season(ts): def load_awi_tracks(): - """All features of the AWI WFS layer (EPSG:3031) with ``km`` added; cached as parquet.""" - pq = CACHE / "awi" / "tracks_3031.parquet" + """AWI ice-sounder features of the WFS layer (EPSG:3031) with ``km``; cached as parquet. + + Fetched in pages of ``AWI_PAGE`` with a server-side filter on + ``radar_system`` (the layer also holds BAS, CReSIS and UTIG data). + """ + pq = CACHE / "awi" / "tracks_3031_sounders.parquet" if pq.exists(): return gpd.read_parquet(pq) - g = gpd.read_file(_fetch(AWI_WFS, CACHE / "awi" / "tracks_3031.geojson")).set_crs(3031, allow_override=True) + cql = "radar_system IN (" + ",".join(f"'{r}'" for r in sorted(AWI_SOUNDERS)) + ")" + pages, start = [], 0 + while True: + q = urllib.parse.urlencode({"service": "WFS", "version": "2.0.0", "request": "GetFeature", + "typeNames": "radar:tracks_3031", "outputFormat": "application/json", + "sortBy": "fid", "count": AWI_PAGE, "startIndex": start, "CQL_FILTER": cql}) + page = gpd.read_file(_fetch(f"{AWI_WFS}?{q}", CACHE / "awi" / f"sounders_{start:06d}.geojson")) + if len(page) == 0: + break + pages.append(page) + start += AWI_PAGE + if len(page) < AWI_PAGE: + break + g = pd.concat(pages, ignore_index=True).set_crs(3031, allow_override=True) g = g[g.geom_type.isin(["LineString", "MultiLineString"])].copy() g["km"] = g.to_crs(4326)["geometry"].apply(geod_km) g.to_parquet(pq) diff --git a/opr_tracks.py b/opr_tracks.py index 53493d5..36e134a 100644 --- a/opr_tracks.py +++ b/opr_tracks.py @@ -10,7 +10,6 @@ """ import re -import urllib.request from pathlib import Path import geopandas as gpd @@ -28,10 +27,8 @@ def _fetch(url, path): - if not path.exists(): - path.parent.mkdir(parents=True, exist_ok=True) - urllib.request.urlretrieve(url, path) - return path + from extra_sources import _fetch as fetch # shared retrying download + return fetch(url, path) def list_segment_csvs(collection): From 1bba79c817daefc4fd12c7ef958e19e0708f534d Mon Sep 17 00:00:00 2001 From: Thomas Teisberg Date: Wed, 16 Sep 2026 14:40:59 -0700 Subject: [PATCH 10/10] Cut density-map memory: per-file BedMap segments (float32) and chunked gridding The CI runner was killed on the density step (~20 GB peak after the shared gap rule loaded all 81M points as one frame). bedmap_segments() now reads one parquet file at a time and bin_line_km() grids in 5M-segment chunks; peak is ~4.8 GB for the full run. Output and the km table are unchanged. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_019UvrC4evHHkT7qWWL96sD5 --- bedmap_campaign_km.csv | 109 +++++++++++++++++++++-------------------- bedmap_common.py | 64 ++++++++++++------------ plot_survey_density.py | 53 +++++++++++--------- 3 files changed, 119 insertions(+), 107 deletions(-) diff --git a/bedmap_campaign_km.csv b/bedmap_campaign_km.csv index bb8108e..a4458b3 100644 --- a/bedmap_campaign_km.csv +++ b/bedmap_campaign_km.csv @@ -1,114 +1,115 @@ name,line_km -BAS_1994_Evans_AIR_BM2,12179.262650741326 -BAS_1998_Dufek_AIR_BM2,4149.363483508933 -BAS_2001_Bailey-Slessor_AIR_BM2,9417.102981384793 +BAS_1994_Evans_AIR_BM2,12179.262650741328 +BAS_1998_Dufek_AIR_BM2,4149.363483508934 +BAS_2001_Bailey-Slessor_AIR_BM2,9417.10298138479 BAS_2001_MAMOG_AIR_BM2,2120.9609719280384 BAS_2001_TORUS_AIR_BM2,7172.650540529199 BAS_2004_BBAS_AIR_BM2,34282.60491095511 -BAS_2005_WISE-ISODYN_AIR_BM2,59688.18212625556 +BAS_2005_WISE-ISODYN_AIR_BM2,59688.18212625557 BAS_2006_GRADES-IMAGE_AIR_BM2,22864.09683745935 BAS_2007_AGAP_AIR_BM2,53245.22499105713 BAS_2007_Lake-Ellsworth_GRN_BM3,291.39961512183146 -BAS_2007_Rutford_GRN_BM3,1274.9219945992415 -BAS_2007_TIGRIS_GRN_BM2,1105.800692939278 +BAS_2007_Rutford_GRN_BM3,1274.9219945992418 +BAS_2007_TIGRIS_GRN_BM2,1105.8006929392782 BAS_2008_Lake-Ellsworth_GRN_BM3,545.0559380028765 -BAS_2009_FERRIGNO_GRN_BM2,1246.2029703357766 +BAS_2009_FERRIGNO_GRN_BM2,1246.2029703357764 BAS_2010_IMAFI_AIR_BM3,21935.345429941688 -BAS_2010_PIG_AIR_BM2,1420.557094373895 +BAS_2010_PIG_AIR_BM2,1420.5570943738949 BAS_2011_Adelaide_AIR_BM3,2496.9405813433536 BAS_2012_Castle_GRN_BM3,212.03544768235724 -BAS_2012_ICEGRAV_AIR_BM3,25010.898751509045 -BAS_2013_ISTAR_GRN_BM3,2558.8971836098713 -BAS_2015_FISS_AIR_BM3,5496.2353004141905 -BAS_2015_POLARGAP_AIR_BM3,28924.060046695704 +BAS_2012_ICEGRAV_AIR_BM3,25010.89875150905 +BAS_2013_ISTAR_GRN_BM3,2558.8971836098717 +BAS_2015_FISS_AIR_BM3,5496.235300414191 +BAS_2015_POLARGAP_AIR_BM3,28924.0600466957 BAS_2016_FISS_AIR_BM3,14099.682974756211 -BAS_2017_English-Coast_AIR_BM3,3719.7530757162654 -BAS_2018_Thwaites_AIR_BM3,4254.490633743226 +BAS_2017_English-Coast_AIR_BM3,3719.753075716266 +BAS_2018_Thwaites_AIR_BM3,4254.490633743227 BAS_2019_Thwaites_AIR_BM3,4496.8047487926015 BEDMAP1_1966-2000_AIR_BM1,512751.4748787665 BGR_1999_GANOVEX-VIII-Matusevich_AIR_BM2,6783.0495726106765 -BGR_1999_GANOVEX-VIII-Mertz_AIR_BM2,1300.4178330820591 +BGR_1999_GANOVEX-VIII-Mertz_AIR_BM2,1300.417833082059 BGR_2002_PCMEGA_AIR_BM2,31769.13425950711 CECS_2006_Subglacial-Lake-CECs_GRN_BM3,1939.2119516739808 -CRESIS_2013_Siple-Coast_AIR_BM3,17274.972127928067 -INGV_1997_ITASE_AIR_BM2,3256.8442284470516 -INGV_1997_Talos-Dome_AIR_BM3,2263.008632785903 -INGV_1999_Talos-Dome_AIR_BM3,11153.019749818553 -INGV_2001_Talos-Dome_AIR_BM3,3817.8735007922646 +CRESIS_2013_Siple-Coast_AIR_BM3,17274.97212792807 +INGV_1997_ITASE_AIR_BM2,3256.844228447052 +INGV_1997_Talos-Dome_AIR_BM3,2263.0086327859035 +INGV_1999_Talos-Dome_AIR_BM3,11153.019749818552 +INGV_2001_Talos-Dome_AIR_BM3,3817.873500792265 INGV_2003_Talos-Dome_AIR_BM3,22061.592422766287 KOPRI_2017_KRT1_AIR_BM3,4339.559925100045 -KOPRI_2018_KRT2_AIR_BM3,1517.4317778335753 +KOPRI_2018_KRT2_AIR_BM3,1517.431777833575 LDEO_2007_AGAP-GAMBIT_AIR_BM2,32609.428173537333 -LDEO_2007_Recovery-Lakes_AIR_BM2,2877.7324600031016 +LDEO_2007_Recovery-Lakes_AIR_BM2,2877.732460003101 LDEO_2015_ROSETTA_AIR_BM3,43523.444895789056 NASA_2002_ICEBRIDGE_AIR_BM2,6238.885501372943 NASA_2004_ICEBRIDGE_AIR_BM2,4290.652034216112 NASA_2009_ICEBRIDGE_AIR_BM2,25425.4062116988 NASA_2010_ICEBRIDGE_AIR_BM2,20071.19778493137 -NASA_2011_ICEBRIDGE_AIR_BM2,65228.757506913884 +NASA_2011_ICEBRIDGE_AIR_BM2,65228.75750691389 NASA_2012_ICEBRIDGE_AIR_BM2,32623.679849093303 -NASA_2013_ICEBRIDGE_AIR_BM3,14166.211692343495 +NASA_2013_ICEBRIDGE_AIR_BM3,14166.211692343497 NASA_2014_ICEBRIDGE_AIR_BM3,52463.17606975831 NASA_2016_ICEBRIDGE_AIR_BM3,55826.74931727524 -NASA_2017_ICEBRIDGE_AIR_BM3,30979.401647501687 -NASA_2018_ICEBRIDGE_AIR_BM3,40577.38595892682 -NASA_2019_ICEBRIDGE_AIR_BM3,20930.313516593418 +NASA_2017_ICEBRIDGE_AIR_BM3,30979.401647501683 +NASA_2018_ICEBRIDGE_AIR_BM3,40577.38595892681 +NASA_2019_ICEBRIDGE_AIR_BM3,20930.313516593415 NIPR_1992_JARE33_GRN_BM3,373.2783262633138 NIPR_1996_JARE37_GRN_BM3,1878.6244437050004 -NIPR_1999_JARE40_GRN_BM3,172.83667006028256 +NIPR_1999_JARE40_GRN_BM3,172.8366700602826 NIPR_2007_JARE49_GRN_BM3,1438.3714194173222 NIPR_2007_JASE_GRN_BM3,2300.734965231171 NIPR_2012_JARE54_GRN_BM3,126.65998551904391 NIPR_2017_JARE59_GRN_BM3,2014.7194283145916 -NIPR_2018_JARE60_GRN_BM3,8905.314937862864 -NPI_2008_BELISSIMA_GRN_BM2,66.30617524699694 -NPI_2010_SRM_AIR_BM2,3270.4030055834273 -NPI_2012_ICERISES_GRN_BM3,327.41267143864115 -NPI_2015_POLARGAP_AIR_BM3,2815.3667153382175 -NPI_2016_MADICE_GRN_BM3,610.904432433296 -PRIC_2004_CHINARE-21_GRN_BM2,2559.982587907926 -PRIC_2007_CHINARE-24_GRN_BM2,603.3034937502166 +NIPR_2018_JARE60_GRN_BM3,8905.314937862866 +NPI_2008_BELISSIMA_GRN_BM2,66.30617524699693 +NPI_2010_SRM_AIR_BM2,3270.403005583427 +NPI_2012_ICERISES_GRN_BM3,327.4126714386412 +NPI_2015_POLARGAP_AIR_BM3,2815.366715338218 +NPI_2016_MADICE_GRN_BM3,610.9044324332962 +PRIC_2004_CHINARE-21_GRN_BM2,2559.9825879079253 +PRIC_2007_CHINARE-24_GRN_BM2,603.3034937502167 PRIC_2015_CHA1_AIR_BM3,40177.03816930688 -PRIC_2016_CHA2_AIR_BM3,21004.283508419878 -PRIC_2017_CHA3_AIR_BM3,25446.214595007812 +PRIC_2016_CHA2_AIR_BM3,21004.28350841988 +PRIC_2017_CHA3_AIR_BM3,25446.21459500781 PRIC_2018_CHA4_AIR_BM3,15052.899018330012 +RNRF_1971_Lambert-Amery_SEI_BM3,0.0 RNRF_1975_Filchner-Ronne_SEI_BM3,1.2835160741144491 -RNRF_1975_Lazarev_SEI_BM3,9.422076160352173 -RNRF_2003_AMSap5_AIR_BM3,4955.185928409161 +RNRF_1975_Lazarev_SEI_BM3,9.422076160352177 +RNRF_2003_AMSap5_AIR_BM3,4955.18592840916 RNRF_2004_AMSap5_AIR_BM3,6360.120969857368 RNRF_2004_Mirny-Vostok_AIR_BM3,1415.5937777432682 -RNRF_2005_AMSap5_AIR_BM3,4247.575813578208 +RNRF_2005_AMSap5_AIR_BM3,4247.575813578207 RNRF_2006_Komsom-Vostok_AIR_BM3,481.9184365624733 RNRF_2006_RAEap5_AIR_BM3,6001.475438810396 -RNRF_2007_AMSap5_AIR_BM3,5729.24149259773 +RNRF_2007_AMSap5_AIR_BM3,5729.241492597731 RNRF_2008_AMSap5_AIR_BM3,5151.903848739583 -RNRF_2009_RAEap5_AIR_BM3,5002.741833256816 -RNRF_2010_RAE_AIR_BM3,7363.233995436161 +RNRF_2009_RAEap5_AIR_BM3,5002.741833256817 +RNRF_2010_RAE_AIR_BM3,7363.233995436162 RNRF_2011_RAE_AIR_BM3,6977.69250118063 -RNRF_2013_RAE_AIR_BM3,14361.90076059565 -RNRF_2013_Vostok-Progress_AIR_BM2,1226.6973135938579 +RNRF_2013_RAE_AIR_BM3,14361.900760595654 +RNRF_2013_Vostok-Progress_AIR_BM2,1226.697313593858 RNRF_2014_RAE_AIR_BM3,6450.733145429845 RNRF_2015_RAE_AIR_BM3,6096.79731173056 RNRF_2016_RAE_AIR_BM3,6490.6332234855 -RNRF_2017_RAE_AIR_BM3,6839.40239536566 +RNRF_2017_RAE_AIR_BM3,6839.402395365659 RNRF_2018_RAE_AIR_BM3,5864.027618573344 RNRF_2019_RAE_AIR_BM3,5654.696345256103 STANFORD_1971_SPRI-NSF-TUD_AIR_BM3,7286.122531415287 STOLAF_1994_Siple-Dome_GRN_BM2,143.5455410337653 STOLAF_2001_ITASE-Byrd-Ellsworth_GRN_BM2,879.0093416878482 STOLAF_2001_ITASE-Ellsworth_GRN_BM2,145.10345224251216 -STOLAF_2002_ITASE-Byrd-South-Pole_GRN_BM2,897.8164427905398 +STOLAF_2002_ITASE-Byrd-South-Pole_GRN_BM2,897.8164427905397 STOLAF_2002_ITASE-Hercules-Dome_GRN_BM2,121.7829483582652 UCANTERBURY_2008_Darwin-Hatherton_GRN_BM2,216.05856338281825 ULB_2012_BEWISE_GRN_BM3,38.48529497266142 ULB_2012_ICECON_GRN_BM3,288.4447429421593 -UTIG_1991_CASERTZ_AIR_BM2,29808.045520122367 -UTIG_1998_West-Marie-Byrd-Land_AIR_BM2,51534.00051425103 +UTIG_1991_CASERTZ_AIR_BM2,29808.04552012237 +UTIG_1998_West-Marie-Byrd-Land_AIR_BM2,51534.00051425102 UTIG_1999_SOAR-LVS-WLK_AIR_BM2,152837.89738200157 UTIG_2000_Robb-Glacier_AIR_BM2,832.2557513557491 -UTIG_2009_Darwin-Hatherton_AIR_BM3,698.5972199718469 -UTIG_2010_ICECAP_AIR_BM3,201552.18361124926 -UTIG_2013_GIMBLE_AIR_BM3,11562.400626045273 +UTIG_2009_Darwin-Hatherton_AIR_BM3,698.5972199718468 +UTIG_2010_ICECAP_AIR_BM3,201552.18361124923 +UTIG_2013_GIMBLE_AIR_BM3,11562.400626045272 UTIG_2015_EAGLE_AIR_BM3,21958.623083374136 -UTIG_2016_OLDICE_AIR_BM3,9142.624259238333 +UTIG_2016_OLDICE_AIR_BM3,9142.624259238331 UWASHINGTON_2018_South-Pole-Lake_GRN_BM3,187.95916516155825 diff --git a/bedmap_common.py b/bedmap_common.py index 66bbb0e..96e0b8e 100644 --- a/bedmap_common.py +++ b/bedmap_common.py @@ -112,48 +112,52 @@ def load_bedmap_catalog(collections=("bedmap2", "bedmap3"), exclude=True): return df.reset_index(drop=True) -def load_bedmap_points(collections=("bedmap1", "bedmap2", "bedmap3"), local_cache=True): - """Kept BedMap points (lon, lat, source_file, row) sorted along track.""" - from xopr.bedmap import fetch_bedmap, query_bedmap - fetch_bedmap() - df = query_bedmap(collections=list(collections), - columns=["lon", "lat", "source_file", "row"], - local_cache=local_cache, show_progress=True) - n_files = df["source_file"].nunique() - df = df[df["source_file"].isin(kept_source_files(collections))] - print(f" {len(df)} BedMap points from {df['source_file'].nunique()} files " - f"({n_files - df['source_file'].nunique()} excluded or duplicate files dropped)") - return df.sort_values(["source_file", "row"]) - - -def point_segments(df, epsg="EPSG:3031", max_point_spacing_m=MAX_POINT_SPACING_M): - """Segments between consecutive points of one file, shorter than the gap limit. - - Returns x1, y1, x2, y2 in ``epsg`` plus the segment lengths in geodesic - km and the source_file of each segment. +def bedmap_segments(epsg="EPSG:3031", collections=("bedmap1", "bedmap2", "bedmap3"), + max_point_spacing_m=MAX_POINT_SPACING_M, with_km=False): + """Gap-filtered flight-line segments of all kept BedMap files. + + Files are fetched to ``radar_cache/bedmap/`` if missing and processed one + at a time (peak memory is one file plus the float32 output, not the 80M + point table). Returns ``(x1, y1, x2, y2)`` in ``epsg`` and, if + ``with_km``, a Series of geodesic km per file. """ + import pyarrow.parquet as pq + import shapely + from xopr.bedmap import fetch_bedmap + paths = sorted(p for v, ps in fetch_bedmap().items() if v in collections for p in ps) + keep = kept_source_files(collections) tf = Transformer.from_crs("EPSG:4326", epsg, always_xy=True) - lon, lat = df["lon"].to_numpy(), df["lat"].to_numpy() - xs, ys = tf.transform(lon, lat) - same_file = df["source_file"].to_numpy()[:-1] == df["source_file"].to_numpy()[1:] - dists = np.hypot(np.diff(xs), np.diff(ys)) - ok = same_file & (dists < max_point_spacing_m) & (dists > 0) - km = _GEOD.inv(lon[:-1][ok], lat[:-1][ok], lon[1:][ok], lat[1:][ok])[2] / 1000.0 - return xs[:-1][ok], ys[:-1][ok], xs[1:][ok], ys[1:][ok], km, df["source_file"].to_numpy()[:-1][ok] + parts, km, n_pts, n_files = [], {}, 0, 0 + for path in paths: + name = Path(path).stem + if name not in keep: + continue + t = pq.read_table(path, columns=["geometry", "row"]).sort_by("row") + c = shapely.get_coordinates(shapely.from_wkb(t.column("geometry").to_numpy(zero_copy_only=False))) + lon, lat = c[:, 0], c[:, 1] + xs, ys = tf.transform(lon, lat) + d = np.hypot(np.diff(xs), np.diff(ys)) + ok = (d < max_point_spacing_m) & (d > 0) + parts.append(np.column_stack([xs[:-1][ok], ys[:-1][ok], xs[1:][ok], ys[1:][ok]]).astype(np.float32)) + if with_km: + km[name] = _GEOD.inv(lon[:-1][ok], lat[:-1][ok], lon[1:][ok], lat[1:][ok])[2].sum() / 1000.0 + n_pts += len(lon); n_files += 1 + seg = np.concatenate(parts) + print(f" {n_pts} BedMap points from {n_files} files ({len(paths) - n_files} excluded or duplicate files dropped)") + return tuple(seg[:, i] for i in range(4)), pd.Series(km, name="line_km").rename_axis("name") def campaign_point_km(rebuild=False): """Gap-filtered line-km per BedMap file, from ``bedmap_campaign_km.csv``. The table is rebuilt from the point data (all three BedMap versions, - ~90M points, needs the point cache) when missing or ``rebuild`` is set; + ~90M points, fetched to the point cache) when missing or ``rebuild`` is set; run ``python bedmap_common.py`` to refresh it after changing the rules. """ if CAMPAIGN_KM_FILE.exists() and not rebuild: return pd.read_csv(CAMPAIGN_KM_FILE, index_col="name")["line_km"] - df = load_bedmap_points() - *_, seg_km, files = point_segments(df) - km = pd.Series(seg_km).groupby(files).sum().rename("line_km").rename_axis("name") + _, km = bedmap_segments(with_km=True) + km = km.sort_index() km.to_csv(CAMPAIGN_KM_FILE) return km diff --git a/plot_survey_density.py b/plot_survey_density.py index 9b997a0..8a15ee7 100644 --- a/plot_survey_density.py +++ b/plot_survey_density.py @@ -71,7 +71,9 @@ def load_bedmap(epsg, local_cache=True, include_extra=True): """BedMap point data (plus extra_sources tracks) as projected segment endpoints. Uses the shared gap rule in ``bedmap_common`` (``MAX_POINT_SPACING_M``), so - the grids and the bar charts count the same line-km. + the grids and the bar charts count the same line-km. Point files are + always read from the local cache (``radar_cache/bedmap/``), fetching any + that are missing; ``local_cache`` is kept for CLI compatibility. Parameters ---------- @@ -86,14 +88,14 @@ def load_bedmap(epsg, local_cache=True, include_extra=True): x1, y1, x2, y2 : np.ndarray Segment endpoint arrays in the target CRS. """ - from bedmap_common import load_bedmap_points, point_segments - df = load_bedmap_points(local_cache=local_cache) - segs = list(point_segments(df, epsg)[:4]) + from bedmap_common import bedmap_segments + segs, _ = bedmap_segments(epsg) + segs = list(segs) if include_extra: from extra_sources import extra_segments ex = extra_segments(epsg) print(f" + {len(ex[0])} segments from extra_sources") - segs = [np.concatenate([a, b]) for a, b in zip(segs, ex)] + segs = [np.concatenate([a, b.astype(np.float32)]) for a, b in zip(segs, ex)] return tuple(segs) def load_xopr(region_filter=None): @@ -136,27 +138,32 @@ def extract_segments(geometries, epsg): np.concatenate(x2s), np.concatenate(y2s)) -def bin_line_km(x1, y1, x2, y2, grid_m, max_extent): - """Bin segment lengths into grid cells, subdividing long segments.""" - dists = np.sqrt((x2 - x1)**2 + (y2 - y1)**2) - ok = dists > 0 - x1, y1, x2, y2, dists = x1[ok], y1[ok], x2[ok], y2[ok], dists[ok] - - n_sub = np.maximum(1, np.ceil(dists / SUBDIV).astype(int)) - total = n_sub.sum() - offsets = np.repeat(np.cumsum(n_sub) - n_sub, n_sub) - idx = np.arange(total) - offsets - rn = np.repeat(n_sub, n_sub).astype(float) - mx = np.repeat(x1, n_sub) + np.repeat(x2 - x1, n_sub) * (idx + 0.5) / rn - my = np.repeat(y1, n_sub) + np.repeat(y2 - y1, n_sub) * (idx + 0.5) / rn - sl = np.repeat(dists, n_sub) / rn +def bin_line_km(x1, y1, x2, y2, grid_m, max_extent, chunk=5_000_000): + """Bin segment lengths into grid cells, subdividing long segments. + Processed in chunks so peak memory stays a few hundred MB for ~80M segments. + """ nx = int(2 * max_extent / grid_m) grid = np.zeros((nx, nx)) - ix = ((mx + max_extent) / grid_m).astype(int) - iy = ((my + max_extent) / grid_m).astype(int) - valid = (ix >= 0) & (ix < nx) & (iy >= 0) & (iy < nx) - np.add.at(grid, (ix[valid], iy[valid]), sl[valid]) + for s in range(0, len(x1), chunk): + a1, b1, a2, b2 = (np.asarray(v[s:s + chunk], dtype=np.float64) for v in (x1, y1, x2, y2)) + dists = np.sqrt((a2 - a1)**2 + (b2 - b1)**2) + ok = dists > 0 + a1, b1, a2, b2, dists = a1[ok], b1[ok], a2[ok], b2[ok], dists[ok] + + n_sub = np.maximum(1, np.ceil(dists / SUBDIV).astype(int)) + total = n_sub.sum() + offsets = np.repeat(np.cumsum(n_sub) - n_sub, n_sub) + idx = np.arange(total) - offsets + rn = np.repeat(n_sub, n_sub).astype(float) + mx = np.repeat(a1, n_sub) + np.repeat(a2 - a1, n_sub) * (idx + 0.5) / rn + my = np.repeat(b1, n_sub) + np.repeat(b2 - b1, n_sub) * (idx + 0.5) / rn + sl = np.repeat(dists, n_sub) / rn + + ix = ((mx + max_extent) / grid_m).astype(int) + iy = ((my + max_extent) / grid_m).astype(int) + valid = (ix >= 0) & (ix < nx) & (iy >= 0) & (iy < nx) + np.add.at(grid, (ix[valid], iy[valid]), sl[valid]) return grid / 1000.0 # m → km