From 884ac82888980bcbfd72de52d2a78de097b296f0 Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Mon, 13 Jul 2026 11:53:48 -0700 Subject: [PATCH 1/4] Add interactive UNR grid web viewer + 0.3 grid data layer Self-contained MapLibre GL viewer (scripts/browse_unr_grid.html) for UNR gridded GPS time series, deployed to GitHub Pages via scripts/deploy-pages.sh (https://opera-adt.github.io/geepers/): - globe (default) / Mercator, 3D terrain, plate boundaries (Bird 2003), theme-aware modern UI (light/dark) - per-point E/N/U time series with +/-sigma band, comparison chart, CSV export; velocity (mm/yr) mode; E+N / Up quiver overlays with scale legend - find-by-id, date-stride subsampling + memory guard, runtime dataset load (file / URL), ETag-keyed fetch to avoid stale cached byte ranges Data pipeline: - create-geoparquet.py: --source grid|stations, --gridded-type, --zero-by mean|start|none, --clear-cache; UTC-aware start_date - UnrGridSource 0.3 (constant/variable gridded products); grid sigma columns converted mm->m like E/N/U - BaseGpsSource: zero_by="none" passthrough; _filter_by_date accepts tz-aware date strings Co-Authored-By: Claude Opus 4.8 --- .gitignore | 9 + README.md | 50 +- scripts/README.md | 121 +- scripts/browse_unr_grid.html | 4716 ++++++++++----------------- scripts/create-geoparquet.py | 184 +- scripts/deploy-pages.sh | 74 + src/geepers/gps_sources/base.py | 142 +- src/geepers/gps_sources/sideshow.py | 20 +- src/geepers/gps_sources/unr.py | 140 +- src/geepers/gps_sources/unr_grid.py | 150 +- src/geepers/utils.py | 23 +- 11 files changed, 2455 insertions(+), 3174 deletions(-) create mode 100755 scripts/deploy-pages.sh diff --git a/.gitignore b/.gitignore index 7163e3a..aacd856 100644 --- a/.gitignore +++ b/.gitignore @@ -161,3 +161,12 @@ Thumbs.db # pixi environments .pixi *.egg-info + +# Local research scratch files (keep out of the package) +/*.zip +/midas.py +/hectorp_wrapper.py +scripts/*.parquet + +# Generated GitHub Pages staging copy (see scripts/deploy-pages.sh) +scripts/web/ diff --git a/README.md b/README.md index 5ad9ee5..564574c 100644 --- a/README.md +++ b/README.md @@ -61,15 +61,55 @@ print(df_many.head()) ``` ## Example: Comparing GPS and InSAR Data -The basic InSAR/GNSS comparison workflow is offered by the `geepers` command line tool. + +The InSAR/GNSS comparison workflow is offered by the `geepers` command line tool: ```bash -geepers --los F33039_los_enu.tif --timeseries-files displacement_20160711_*tif --temporal-coherence-files temporal_coherence_*.tif --similarity-files phase_similarity*.tif +geepers \ + --los F33039_los_enu.tif \ + --timeseries-files displacement_20160711_*tif \ + --temporal-coherence-files temporal_coherence_*.tif \ + --similarity-files phase_similarity*.tif \ + --insar-buffer-meters 100 \ + --requirement-mm 3 0.5 \ + --wavelength 0.2384 # NISAR L-band; omit for Sentinel-1 ``` -The results are saved in the current directory in the `GPS` folder by default. - -(TODO: Example data prep for this) +Results are saved in the `GPS` folder by default: per-station time-series +and rate comparisons, plus the structure function (pairwise relative RMSE +vs station separation, checked against the requirement curve) and the +per-epoch network misfit. See +[How-To Guides](https://geepers.readthedocs.io/en/latest/how-to-guides/) +and the runnable +[validation notebook](docs/notebooks/gnss_insar_validation.ipynb). + +## Analysis toolbox + +Beyond data access and the comparison workflow, geepers ships analysis +modules for GNSS velocity fields (see +[Analysis modules](docs/analysis-modules.md) for usage and the +[tour notebook](docs/notebooks/geepers_tour.ipynb) for a runnable demo): + +| Module | Purpose | +|---|---| +| `geepers.midas` | Robust MIDAS velocities (Blewitt et al., 2016) | +| `geepers.trend` | Velocities with realistic uncertainties under power-law + white noise (validated against HectorP); fast Whittle method and parallel `estimate_trend_many` for networks | +| `geepers.variability` | Temporal & spatial velocity-stability metrics, spatial structure function | +| `geepers.quality` | Gap percentage, station quality, reference selection | +| `geepers.steps` | Detection of uncatalogued jumps (AIC sliding window) | +| `geepers.cme` | Common-mode error estimation/removal (PCA/ICA) | +| `geepers.gps_imaging` | Robust weighted-median interpolation (Hammond et al., 2016 GPS Imaging port) | +| `geepers.collocation` | Least-squares collocation, ordinary kriging, plate-boundary separation | +| `geepers.euler` | Euler pole estimation and plate-motion prediction | +| `geepers.strain` | Strain-rate/rotation fields from gridded velocities | +| `geepers.validation` | GNSS-vs-InSAR validation: velocity scatter (MAD/RMSE/R²), structure function, semivariogram, per-epoch misfit | +| `geepers.synthetic` | Schema-valid synthetic networks and series for testing | + +An interactive MapLibre viewer for UNR gridded time series is hosted at +**[opera-adt.github.io/geepers](https://opera-adt.github.io/geepers/)** +(globe view, velocity/vector overlays, plate boundaries, per-point time +series). See the [viewer docs](scripts/README.md) to run it locally or +build your own dataset. ### Working with Multiple Sources diff --git a/scripts/README.md b/scripts/README.md index d8ea4c9..e01dfa3 100644 --- a/scripts/README.md +++ b/scripts/README.md @@ -1,13 +1,120 @@ -# UNR Grid Web Browser +# OPERA UNR Grid Web Browser -Setup: +Interactive MapLibre GL viewer for UNR gridded GPS time series. +The gridded data are produced by the Nevada Geodetic Laboratory (UNR), +funded by the JPL-led [OPERA](https://www.jpl.nasa.gov/go/opera) project; +the viewer is developed at JPL. + +> **Disclaimer**: the viewer and the underlying gridded GPS products are +> research tools provided "as is", without warranty of any kind. +> Displacements, uncertainties, and derived velocities are experimental +> and may contain errors or artifacts. Use of this tool does not imply +> endorsement by JPL/Caltech, NASA, or the University of Nevada, Reno. +Loads a single Parquet file directly in the browser (via +[hyparquet](https://github.com/hyparam/hyparquet)) and scrubs through dates +with GPU-driven color updates — no per-date files, no re-fetching. + +## Live site + +Hosted on GitHub Pages: **https://opera-adt.github.io/geepers/** + +The default dataset is the **global UNR grid, all 28,358 points at monthly +sampling** (2014→2026, ~95 MB), served same-origin from the `gh-pages` +branch. See `deploy-pages.sh` for how the site is (re)built and pushed. + +### Viewing the full daily-resolution grid + +The full daily grid is too large for browser hosting — GitHub has no +surface that serves it cross-origin (Pages caps files at 100 MB; Release +assets and LFS send no CORS headers). It ships instead as a **local +artifact** for offline viewing: ```bash cd scripts/ -# EXAMPLE BBOX: -python create-geojson.py --bbox -110 28 -101 36 --start-date 2016-01-01 -# Creates geojson_sources/ -echo "Visit the URL in your browser:" -echo "http://localhost:8123/browse_unr_grid.html" +# OPERA_UNR_GNSS_grid_full.parquet: all 28,358 points, daily, 2014→2026 +# (~850 MB; viewer-minimal columns date_idx/point_idx/E/N/U, no sigmas) python -m http.server 8123 +# Open http://localhost:8123/browse_unr_grid.html and use "Open .parquet…" +# in the Data panel to pick OPERA_UNR_GNSS_grid_full.parquet, or: +# browse_unr_grid.html?data=OPERA_UNR_GNSS_grid_full.parquet ``` + +At full daily resolution the viewer needs ~1.5 GB of browser memory (it +will ask to confirm); use the **Date stride** selector (or `?stride=N`) to +subsample and lighten it. The `_full` file omits the sigma columns, so the +±σ chart band is unavailable there — use a smaller/regional export (with +sigmas) if you need uncertainties. + +## Setup + +```bash +cd scripts/ +# 1. Download data and build the viewer-ready Parquet file (example bbox): +python create-geoparquet.py --bbox -110 28 -101 36 --start-date 2016-01-01 +# Creates unr_grid.parquet + +# 2. Serve and open: +python -m http.server 8123 +# Visit http://localhost:8123/browse_unr_grid.html +``` + +Useful options: + +- `--source grid|stations` — UNR gridded (interpolated) product, or real + UNR GPS station positions (.tenv3). Default `grid`. +- `--gridded-type constant|variable` — time-constant vs time-variable UNR + product (version 0.3 only; grid source only; default `variable`). +- `--output-file my_area.parquet` then open + `browse_unr_grid.html?data=my_area.parquet`. +- `--clear-cache` — wipe the geepers download cache first (forces fresh + downloads). +- `--zero-by mean|start|none` — zero each point's series by its mean, its + first epochs, or `none` to keep values exactly as published (default + `mean`). +- If no file is found, the page offers a local file picker (drag any + compatible `.parquet` in — nothing is uploaded, parsing is in-browser). + +The viewer is a single self-contained HTML file — MapLibre GL v5, uPlot and +hyparquet are inlined, so only the basemap/terrain tiles need the network. + +## Viewer features + +- Date slider + playback (2–30 fps), keyboard: `←`/`→` step, `space` play. +- Click a grid point → East/North/Up time series chart (uPlot) with an + optional ±1σ shaded band; **Shift+click** a second point for a comparison + chart. Charts are resizable (drag the corner), zoomable (drag box, mouse + wheel, double-click resets), and clicking a sample jumps the map to that + date. +- Component selector, colormaps (RdBu, BrBG, Viridis, Turbo, Magma), + invert, symmetric/robust (p2–p98) or manual range in mm; `live` re-runs + the auto range on every date change while scrubbing/playing. +- Velocity mode: color points by per-point linear trend (least-squares, + mm/yr) instead of per-date displacement. +- Vector overlay: horizontal (E+N) and/or vertical (Up, red up / blue + down) quiver arrows over the points, with an arrow-scale slider and a + scale legend above the colorbar. Arrow scaling is automatic (p90 of the + data, follows the date like the color `live` mode) or fixed via typed + reference magnitudes (e.g. H 3, V 1 mm/yr). Arrows show the same field + as the colors (per-date displacement, or velocity in velocity mode). + The globe view gets a dark space backdrop. +- Each chart has a `csv` button (dates + E/N/U ± σ of that point, in mm, + with the current referencing applied). +- Find ID box zooms to a grid point / station by identifier. +- Large files: a memory estimate is checked before loading, and a "Date + stride" selector (`?stride=N`) loads only every Nth date to bound + memory (e.g. the ~1 GB time-variable CA file fits comfortably with + stride 5). +- Reference modes: none (values exactly as stored in the file), per-point + temporal mean, first date, or any chosen date (displacement relative to + that date). +- Basemaps: Carto light/dark, OSM, Esri satellite; globe (default) or + Mercator projection; optional 3D terrain (AWS terrain tiles) with + adjustable exaggeration — right-drag / Ctrl+drag to tilt and rotate. +- Tectonic plate boundaries overlay (Bird 2003, via + [fraxen/tectonicplates](https://github.com/fraxen/tectonicplates)); + loads `PB2002_boundaries.json` next to the HTML if present, else from + GitHub raw. +- Data panel: load another `.parquet` (local file or URL) without reloading + the page, and a "Clear cache & reload" button. Fetches are keyed to the + file's `Last-Modified`/`ETag`, so regenerating a parquet under the same + name can never serve stale cached byte ranges. diff --git a/scripts/browse_unr_grid.html b/scripts/browse_unr_grid.html index bc08708..7e4ea31 100644 --- a/scripts/browse_unr_grid.html +++ b/scripts/browse_unr_grid.html @@ -4,3194 +4,1924 @@ - UNR Gridded Time Series Viewer + OPERA UNR Gridded Time Series Viewer + - - - - - - - - + + + + - - - - - - - -
+ - -
- -
- - - +
+ +
Loading…
+
+
+ +
+

OPERA UNR Grid Viewer

+
+ Gridded GPS time series by the Nevada Geodetic Laboratory (UNR), + funded by the JPL-led OPERA project. Viewer: JPL.
- -
- -
-
-
- - Map Controls -
-
- - -
- - - -
- - - - -
- -
- - -
-
- - -
- - -
- - -
- - - -
- -
-
- - -
- - -
- - -
-
- - -
-
-
- - -
- -
-
- - -
-
- - -
-
-
-
-
Visible: -
-
Valid: -
-
Min: -
-
Max: -
-
-
- Statistics for current view extent -
-
-
-
+
+ Display +
+ +
- - -
-
-
- - Performance Settings -
-
- -
- -
- -
- Max cache size: 50 files - -
-
-
- -
- -
- -
-
- -
- -
- - -
- Controls data points loaded for popup charts. "All data" may be slow for long time ranges. -
-
-
- -
- -
- -
- Displays cache status, load times, and render times -
-
-
- -
-
-
- - Performance Tips -
-
    -
  • • Enable caching for smoother playback
    • -
  • • Use viewport culling for large datasets
    • -
  • • Lower cache size if memory is limited
    • -
  • • Use "Smart sampling" for time series unless you need all data points
    • -
  • • "All data" mode provides complete time series but loads slower
    • -
  • • Monitor metrics to optimize settings
    • -
-
-
+
+ +
- - -
-
-
- - Map Animation Export -
-
- -
- -
-
- - -
-
- - -
-
-
- - -
-
- -
- -
-
- - -
-
- - -
-
-
- - -
-
- -
- -
- - - - -
-
- -
- -
-
-
Preview:
-
-
Dates: -
-
Frames: -
-
Duration: -
-
File size: ~-
-
-
- - - - - - - - -
-
- -
-
-
- - How Map Animation Export Works -
-
-
Rendering: Creates a basic map representation with data points
-
High Quality: Records at chosen resolution with smooth animation
-
Customizable: Control date range, speed, and visual elements
-
Professional: Creates smooth MP4/WebM videos ready for presentations
-
-
- - Note: Export captures a simplified map representation rather than full tile imagery due to browser security limitations. -
-
-
+
+ + +
-
- - - - - - - - - + + + + - // Show all hover elements - infoBg.style.opacity = '1'; - infoDate.style.opacity = '1'; - infoValue.style.opacity = '1'; - hoverLine.style.opacity = '0.7'; + - // Highlight the specific point - dataPoints.forEach(dp => { - dp.setAttribute('r', '3'); - dp.setAttribute('fill', '#3b82f6'); - }); - dataPoints[index].setAttribute('r', '5'); - dataPoints[index].setAttribute('fill', '#1d4ed8'); - }; - - const hideInfo = () => { - infoBg.style.opacity = '0'; - infoDate.style.opacity = '0'; - infoValue.style.opacity = '0'; - hoverLine.style.opacity = '0'; - - // Reset all points - dataPoints.forEach(dp => { - dp.setAttribute('r', '3'); - dp.setAttribute('fill', '#3b82f6'); - }); - }; + diff --git a/scripts/create-geoparquet.py b/scripts/create-geoparquet.py index 299d77d..1e2cfc1 100644 --- a/scripts/create-geoparquet.py +++ b/scripts/create-geoparquet.py @@ -1,72 +1,190 @@ +"""Export UNR gridded time series to a browser-optimized Parquet file. + +The output file is consumed by `browse_unr_grid.html` (MapLibre GL viewer), +and is equally usable from pandas / GeoPandas / DuckDB: + + duckdb -c "SELECT * FROM 'unr_grid.parquet' LIMIT 5" + +Layout notes +------------ +- Long format: one row per (grid point, date). +- Sorted by (date, point) and written with snappy compression, which + hyparquet can decompress natively in the browser (no extra codecs). +- `date_idx` / `point_idx` integer columns let the viewer scatter values + into dense [n_dates x n_points] matrices without parsing dates or ids. +- The full date list and point coordinates are embedded as JSON in the + Parquet file-level metadata (key ``unr_grid_meta``), so the viewer can + build the map without scanning string columns. +""" + import datetime import json +import shutil from pathlib import Path from typing import Literal +import numpy as np import pandas as pd +import pyarrow as pa +import pyarrow.parquet as pq import tyro -from geepers.gps_sources import UnrGridSource +from geepers.gps_sources import BaseGpsSource, UnrGridSource, UnrSource +META_KEY = "unr_grid_meta" -def export_gdf_to_geoparquet(gdf, output_file="unr_grid.parquet"): - """Export a GeoDataFrame to multiple GeoJSON files organized by date. + +def export_gdf_to_parquet(gdf, output_file="unr_grid.parquet") -> Path: + """Export a long-format GeoDataFrame from `timeseries_many` to Parquet. Parameters ---------- gdf : GeoDataFrame - Result from `timeseries_many` - output_dir : str - Directory to save GeoJSON files + Result from `timeseries_many` (one row per point per date). + output_file : str | Path + Output .parquet path. Returns ------- - dict - Mapping of source names to file paths + Path + Path to the written file. """ - # Ensure date column is datetime - if not pd.api.types.is_datetime64_any_dtype(gdf["date"]): - gdf["date"] = pd.to_datetime(gdf["date"]) - - # Create output directory - output_path = Path(output_file) - if output_path.suffix.lower() != ".parquet": - output_path = output_path.with_suffix(".parquet") + output_path = Path(output_file).with_suffix(".parquet") output_path.parent.mkdir(parents=True, exist_ok=True) - gdf.to_parquet( - output_file, compression="snappy", row_group_size=None, geometry_encoding="WKB" + df = pd.DataFrame(gdf.drop(columns="geometry", errors="ignore")) + if not pd.api.types.is_datetime64_any_dtype(df["date"]): + df["date"] = pd.to_datetime(df["date"]) + + df = df.sort_values(["date", "id"], kind="mergesort", ignore_index=True) + + # Integer indices for fast dense-matrix assembly in the browser + dates = df["date"].dt.normalize() + unique_dates = dates.drop_duplicates().reset_index(drop=True) + df["date_idx"] = dates.map( + pd.Series(np.arange(len(unique_dates), dtype=np.int32), index=unique_dates) + ) + point_codes, unique_ids = pd.factorize(df["id"], sort=True) + df["point_idx"] = point_codes.astype(np.int32) + + points = ( + df.drop_duplicates("point_idx") + .sort_values("point_idx")[["id", "lon", "lat"]] + .reset_index(drop=True) + ) + + value_cols = ["east", "north", "up", "sigma_east", "sigma_north", "sigma_up"] + out = pd.DataFrame( + { + "id": df["id"].astype("string"), + "date": df["date"].dt.date, # date32: 4 bytes, no tz ambiguity + "date_idx": df["date_idx"], + "point_idx": df["point_idx"], + "lon": df["lon"].astype(np.float32), + "lat": df["lat"].astype(np.float32), + **{c: df[c].astype(np.float32) for c in value_cols}, + } + ) + + meta = { + "dates": [d.strftime("%Y-%m-%d") for d in unique_dates], + "points": { + "id": points["id"].tolist(), + "lon": [round(float(v), 6) for v in points["lon"]], + "lat": [round(float(v), 6) for v in points["lat"]], + }, + "value_columns": value_cols, + "units": "meters", + } + + table = pa.Table.from_pandas(out, preserve_index=False) + table = table.replace_schema_metadata( + {**(table.schema.metadata or {}), META_KEY.encode(): json.dumps(meta).encode()} + ) + # Row groups aligned to whole dates keep per-date reads contiguous + n_points = len(points) + rows_per_group = max(n_points * max(1, 262_144 // max(n_points, 1)), n_points) + pq.write_table( + table, + output_path, + compression="snappy", # hyparquet decodes snappy without extra codecs + row_group_size=rows_per_group, + use_dictionary=["id"], ) - # Print file size statistics - file_size_mb = output_path.stat().st_size / (1024 * 1024) - print(f"\nOutput file size: {file_size_mb:.2f} MB") - print(f"Successfully created: {output_file}") + size_mb = output_path.stat().st_size / 2**20 + print( + f"Wrote {output_path} ({size_mb:.1f} MB): " + f"{len(out):,} rows, {n_points} points, {len(unique_dates)} dates" + ) + return output_path def main( bbox: tuple[float, float, float, float], + source: Literal["grid", "stations"] = "grid", start_date: datetime.datetime = datetime.datetime(2016, 1, 1), - output_dir=Path("geojson_sources"), - version: Literal["0.1", "0.2"] = "0.2", + output_file: Path = Path("unr_grid.parquet"), + version: Literal["0.1", "0.3"] = "0.3", + gridded_type: Literal["constant", "variable"] = "variable", + cache_dir: Path | None = None, + max_workers: int = 8, + clear_cache: bool = False, + zero_by: Literal["mean", "start", "none"] = "mean", ): - """Export a GeoDataFrame to multiple GeoJSON files organized by date. + """Download UNR time series and export a viewer-ready Parquet file. Parameters ---------- bbox : tuple[float, float, float, float] - Result from `timeseries_many` - output_dir : str - Directory to save GeoJSON files + Bounding box (west, south, east, north) in degrees. + source : {"grid", "stations"} + "grid" downloads the UNR gridded (interpolated) product; + "stations" downloads real UNR GPS station positions (.tenv3). + Default is "grid". start_date : datetime - First date to download from UNR. - Default is 2016-01-01 + First date to keep. Default is 2016-01-01. + output_file : Path + Output .parquet path. Default is unr_grid.parquet. + version : {"0.1", "0.3"} + UNR grid data version (grid source only). + gridded_type : {"constant", "variable"} + Time-constant or time-variable gridded product (0.3 only; grid source only). + cache_dir : Path, optional + Where downloaded .tenv8/.tenv3 files are cached. + Default is ~/.cache/geepers. + max_workers : int + Parallel download threads. Default is 8. + clear_cache : bool + Delete this source's download cache before fetching, forcing + fresh downloads. Default is False. + zero_by : {"mean", "start", "none"} + How each point's time series is zeroed: subtract its mean, its + first ~10 epochs, or "none" to keep values exactly as published. + Default is "mean". """ - unrg = UnrGridSource(version=version) - gdf = unrg.timeseries_many(bbox=bbox, start_date=start_date) - export_gdf_to_geoparquet(gdf=gdf, output_dir=output_dir) + src: BaseGpsSource + if source == "grid": + src = UnrGridSource( + version=version, gridded_type=gridded_type, cache_dir=cache_dir + ) + else: + src = UnrSource(cache_dir=cache_dir) + + if clear_cache: + print(f"Clearing download cache: {src._cache_dir}") + shutil.rmtree(src._cache_dir, ignore_errors=True) + src._cache_dir.mkdir(parents=True, exist_ok=True) + + gdf = src.timeseries_many( + bbox=bbox, + start_date=start_date.replace(tzinfo=datetime.UTC).isoformat(), + zero_by=zero_by, + max_workers=max_workers, + ) + export_gdf_to_parquet(gdf=gdf, output_file=output_file) if __name__ == "__main__": diff --git a/scripts/deploy-pages.sh b/scripts/deploy-pages.sh new file mode 100755 index 0000000..03d29a8 --- /dev/null +++ b/scripts/deploy-pages.sh @@ -0,0 +1,74 @@ +#!/usr/bin/env bash +# Rebuild and deploy the GitHub Pages site for the OPERA UNR grid viewer. +# +# Publishes a single self-contained page plus its default dataset to the +# `gh-pages` branch of opera-adt/geepers, served at +# https://opera-adt.github.io/geepers/ +# +# The branch is rebuilt as ONE fresh orphan commit each run so old (large) +# data blobs never accumulate in history. +# +# Usage: +# ./deploy-pages.sh [DATA_PARQUET] +# DATA_PARQUET defaults to OPERA_UNR_GNSS_grid_monthly.parquet (the global +# monthly grid). It must be < 100 MB (GitHub Pages per-file limit). +set -euo pipefail + +REPO="opera-adt/geepers" +REMOTE="https://github.com/${REPO}" +SCRIPTS_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" + +DATA_SRC="${1:-${SCRIPTS_DIR}/OPERA_UNR_GNSS_grid_monthly.parquet}" +# The .zip extension is deliberate: it stops the GitHub Pages CDN from +# gzipping the file, which would corrupt hyparquet's HTTP range reads. +DATA_DEST="OPERA_UNR_GNSS_grid.parquet.zip" +BOUNDARIES="${SCRIPTS_DIR}/PB2002_boundaries.json" + +[ -f "$DATA_SRC" ] || { echo "error: data file not found: $DATA_SRC" >&2; exit 1; } +size_mb=$(( $(stat -c%s "$DATA_SRC") / 1024 / 1024 )) +if [ "$size_mb" -ge 100 ]; then + echo "error: $DATA_SRC is ${size_mb} MB; GitHub Pages caps files at 100 MB." >&2 + exit 1 +fi + +WORK="$(mktemp -d)" +trap 'rm -rf "$WORK"' EXIT + +# 1. Viewer: point the default DATA_URL at the hosted (renamed) dataset. +python3 - "$SCRIPTS_DIR/browse_unr_grid.html" "$WORK/index.html" "$DATA_DEST" <<'PY' +import sys +src, dst, data = sys.argv[1:4] +html = open(src).read() +old = "const DATA_URL = params.get('data') || 'unr_grid.parquet';" +new = f"const DATA_URL = params.get('data') || '{data}';" +assert html.count(old) == 1, "could not find DATA_URL default in viewer" +open(dst, "w").write(html.replace(old, new)) +PY + +# 2. Static assets. +cp "$DATA_SRC" "$WORK/$DATA_DEST" +[ -f "$BOUNDARIES" ] && cp "$BOUNDARIES" "$WORK/PB2002_boundaries.json" +: > "$WORK/.nojekyll" +cat > "$WORK/README.md" <\` (host must allow CORS + ranges). +EOF + +# 3. Fresh single-commit orphan branch, force-pushed. +cd "$WORK" +git init -q -b gh-pages +git add -A +git commit -q -m "Deploy OPERA UNR grid viewer ($(basename "$DATA_SRC"), ${size_mb} MB)" +git remote add origin "$REMOTE" +git push -f origin gh-pages + +echo "Deployed. Pages will rebuild shortly: https://opera-adt.github.io/geepers/" diff --git a/src/geepers/gps_sources/base.py b/src/geepers/gps_sources/base.py index cdb6511..ea532ef 100644 --- a/src/geepers/gps_sources/base.py +++ b/src/geepers/gps_sources/base.py @@ -3,6 +3,9 @@ from __future__ import annotations import difflib +import logging +import re +import urllib.error import warnings from abc import ABC, abstractmethod from collections.abc import Iterable @@ -11,7 +14,7 @@ import geopandas as gpd import pandas as pd -from shapely.geometry import box +import requests from tqdm.auto import tqdm from tqdm.contrib.concurrent import thread_map @@ -19,7 +22,39 @@ from geepers._types import PathOrStr from geepers.schemas import PointSchema -__all__ = ["BaseGpsSource"] +__all__ = ["BaseGpsSource", "validate_station_id"] + +logger = logging.getLogger("geepers") + +_STATION_ID_PATTERN = re.compile(r"[A-Za-z0-9_]{1,9}") + + +def validate_station_id(station_id: str) -> str: + """Check that `station_id` is safe to embed in URLs and cache paths. + + Guards against path traversal (e.g. ``"../../etc"``), since station ids + are interpolated into both download URLs and local cache filenames. + + Parameters + ---------- + station_id : str + The station identifier to check. + + Returns + ------- + str + The validated station id, unchanged. + + Raises + ------ + ValueError + If the id contains anything but 1-9 alphanumeric/underscore chars. + + """ + if not _STATION_ID_PATTERN.fullmatch(station_id): + msg = f"Invalid station id: {station_id!r}" + raise ValueError(msg) + return station_id class BaseGpsSource(ABC): @@ -34,39 +69,43 @@ def timeseries_many( frame: Literal["ENU", "XYZ"] = "ENU", start_date: str | None = None, end_date: str | None = None, - zero_by: Literal["mean", "start"] = "mean", + zero_by: Literal["mean", "start", "none"] = "mean", download_if_missing: bool = True, *, max_workers: int = 8, + skip_errors: bool = True, ): if bbox is None and mask is None and ids is None: msg = "Must provide ids, bbox or mask" raise ValueError(msg) - gdf_stations = self.stations() - if bbox is not None: - gdf_stations = self.stations(bbox=bbox) - ids = gdf_stations["id"] - elif mask is not None: - gdf_stations = self.stations(mask=mask) - ids = gdf_stations["id"] - - if ids is None: + gdf_stations = self.stations(bbox=bbox, mask=mask) + if bbox is not None or mask is not None or ids is None: ids = gdf_stations["id"] + # Index by id for O(1) metadata lookups in `_load_one` + station_rows = gdf_stations.set_index("id") # Function to load one id - def _load_one(sid: str) -> pd.DataFrame: - df = self.timeseries( - sid, - frame=frame, - start_date=start_date, - end_date=end_date, - zero_by=zero_by, - download_if_missing=download_if_missing, - ) + def _load_one(sid: str) -> pd.DataFrame | None: + try: + df = self.timeseries( + sid, + frame=frame, + start_date=start_date, + end_date=end_date, + zero_by=zero_by, + download_if_missing=download_if_missing, + ) + except (requests.HTTPError, urllib.error.HTTPError) as e: + # Some ids in the station lists have no data file on the + # server (e.g. UNR grid points with too little data) + if not skip_errors: + raise + logger.warning("Skipping %s: %s", sid, e) + return None df.insert(0, "id", sid) # keep id as a column for melt/pivot - row = gdf_stations[gdf_stations["id"] == sid] + row = station_rows.loc[sid] for col in ("lon", "lat", "alt", "geometry"): - df[col] = row.iloc[0][col] + df[col] = row[col] return df # (Optional) parallel map @@ -77,6 +116,13 @@ def _load_one(sid: str) -> pd.DataFrame: else: dfs = [_load_one(sid) for sid in tqdm(ids)] + n_failed = sum(df is None for df in dfs) + if n_failed: + logger.warning("Skipped %d of %d ids with no data", n_failed, len(dfs)) + dfs = [df for df in dfs if df is not None] + if not dfs: + msg = "No time series could be loaded" + raise ValueError(msg) big = pd.concat(dfs, ignore_index=True) return gpd.GeoDataFrame(big, geometry="geometry", crs="EPSG:4326") @@ -96,9 +142,18 @@ def __init__(self, cache_dir: PathOrStr | None = None): self._base_cache_dir = utils.get_cache_dir() else: self._base_cache_dir = Path(cache_dir) - subdir = self.__class__.__name__.lower().replace("source", "") - self._cache_dir = self._base_cache_dir / subdir - self._cache_dir.mkdir(exist_ok=True, parents=True) + self._subdir = self.__class__.__name__.lower().replace("source", "") + + @property + def _cache_dir(self) -> Path: + """Cache directory for this source, created on first access. + + Creation is deferred so that instantiating a source (e.g. at module + import time) has no filesystem side effects. + """ + cache_dir = self._base_cache_dir / self._subdir + cache_dir.mkdir(exist_ok=True, parents=True) + return cache_dir @abstractmethod def timeseries( @@ -108,7 +163,7 @@ def timeseries( frame: Literal["ENU", "XYZ"] = "ENU", start_date: str | None = None, end_date: str | None = None, - zero_by: Literal["mean", "start"] = "mean", + zero_by: Literal["mean", "start", "none"] = "mean", download_if_missing: bool = True, ) -> pd.DataFrame: """Load GPS station time series data. @@ -197,18 +252,18 @@ def _apply_spatial_filters( Filtered GeoDataFrame. """ - # Apply bbox filter + # Apply bbox filter (coordinate slicing: much faster than a + # geometric clip for point layers) if bbox is not None: west, south, east, north = bbox - bounds_poly = box(west, south, east, north) - gdf = gdf.clip(bounds_poly) + gdf = gdf.cx[west:east, south:north] # Apply mask filter if mask is not None: - gdf = gdf[gdf.geometry.within(mask.unary_union)] + gdf = gdf[gdf.geometry.within(mask.union_all())] # Reset index for cleaner output - gdf.reset_index(drop=True, inplace=True) + gdf = gdf.reset_index(drop=True) # Validate basic point schema return PointSchema.validate(gdf, lazy=True) @@ -219,22 +274,33 @@ def _filter_by_date( start_date: str | None = None, end_date: str | None = None, ) -> pd.DataFrame: - """Filter DataFrame by date range.""" + """Filter DataFrame by date range. + + Accepts tz-aware date strings (e.g. UTC ISO timestamps); the tz is + dropped so the bound compares against the tz-naive ``date`` column. + """ + + def _naive(value: str) -> pd.Timestamp: + ts = pd.to_datetime(value) + return ts.tz_localize(None) if ts.tzinfo is not None else ts + if start_date: - df = df[df["date"] >= pd.to_datetime(start_date)] + df = df[df["date"] >= _naive(start_date)] if end_date: - df = df[df["date"] <= pd.to_datetime(end_date)] + df = df[df["date"] <= _naive(end_date)] return df def _zero_data( self, df: pd.DataFrame, - zero_by: Literal["mean", "start"] = "mean", + zero_by: Literal["mean", "start", "none"] | None = "mean", columns: list[str] | None = None, ) -> pd.DataFrame: - """Zero the data in a DataFrame.""" + """Zero the data in a DataFrame ("none"/None leaves it as-is).""" if columns is None: columns = ["east", "north", "up"] + if zero_by is None or zero_by.lower() == "none": + return df if zero_by.lower() == "mean": mean_val = df[columns].mean() df.loc[:, columns] -= mean_val @@ -242,7 +308,7 @@ def _zero_data( start_val = df[columns].iloc[:10].mean() df.loc[:, columns] -= start_val else: - msg = "zero_by must be either 'mean' or 'start'" + msg = "zero_by must be 'mean', 'start', or 'none'" raise ValueError(msg) return df diff --git a/src/geepers/gps_sources/sideshow.py b/src/geepers/gps_sources/sideshow.py index 560129d..4ff8fdc 100644 --- a/src/geepers/gps_sources/sideshow.py +++ b/src/geepers/gps_sources/sideshow.py @@ -3,18 +3,15 @@ from __future__ import annotations import logging -import warnings from functools import lru_cache -from io import StringIO from typing import TYPE_CHECKING, Final, Literal import geopandas as gpd import numpy as np import pandas as pd -import requests from geepers.schemas import StationObservationSchema -from geepers.utils import decimal_year_to_datetime, get_cache_dir +from geepers.utils import decimal_years_to_datetimes, get_cache_dir from .base import BaseGpsSource @@ -62,7 +59,7 @@ def timeseries( frame: Literal["ENU", "XYZ"] = "ENU", start_date: str | None = None, end_date: str | None = None, - zero_by: Literal["mean", "start"] = "mean", + zero_by: Literal["mean", "start", "none"] = "mean", download_if_missing: bool = True, # noqa: ARG002 ) -> pd.DataFrame: """Load GPS station time series data. @@ -92,9 +89,10 @@ def timeseries( msg = "XYZ frame not supported for Sideshow data" raise ValueError(msg) - df = self._read_series(station_id) + # Copy so we never mutate the lru_cached DataFrame + df = self._read_series(station_id).copy() # Replace decimal year with datetime - df["date"] = df["decimal_year"].apply(decimal_year_to_datetime) + df["date"] = decimal_years_to_datetimes(df["decimal_year"]) df = df.drop(columns=["decimal_year"]) # Move date to first column: df = df[["date", *df.columns[:-1].to_list()]] @@ -103,7 +101,7 @@ def timeseries( return StationObservationSchema.validate(df, lazy=True) @staticmethod - @lru_cache(maxsize=1) + @lru_cache(maxsize=128) def _read_series(station_id: str) -> pd.DataFrame: _raw_names = ["decimal_year"] + SideshowSource._names[1:] # https://sideshow.jpl.nasa.gov/post/tables/GNSS_Time_Series.pdf @@ -128,10 +126,10 @@ def _read_series(station_id: str) -> pd.DataFrame: @lru_cache(maxsize=1) def _fetch_station_data() -> gpd.GeoDataFrame: """Download and cache the JPL Sideshow site list.""" - resp = requests.get(SITE_LIST_URL) - resp.raise_for_status() + import warnings + with warnings.catch_warnings(category=UserWarning, action="ignore"): - return np.loadtxt(StringIO(resp.text), comments="<", skiprows=9, dtype=str) + return np.loadtxt(SITE_LIST_URL, comments="<", skiprows=9, dtype=str) def _read_station_data(self) -> gpd.GeoDataFrame: lines = self._fetch_station_data() diff --git a/src/geepers/gps_sources/unr.py b/src/geepers/gps_sources/unr.py index 1545e28..25ad40b 100644 --- a/src/geepers/gps_sources/unr.py +++ b/src/geepers/gps_sources/unr.py @@ -5,7 +5,6 @@ import datetime import logging from functools import cache -from pathlib import Path from typing import Literal import geopandas as gpd @@ -13,11 +12,10 @@ import pandas as pd import requests -from geepers import utils from geepers._types import PathOrStr from geepers.schemas import StationObservationSchema -from .base import BaseGpsSource +from .base import BaseGpsSource, validate_station_id __all__ = ["UnrSource"] @@ -30,10 +28,11 @@ # https://geodesy.unr.edu/gps_timeseries/IGS20/tenv3/IGS20/LAVR.tenv3 "https://geodesy.unr.edu/gps_timeseries/{reference}/tenv3/{reference}/{station}.tenv3" ) -GPS_DIR = utils.get_cache_dir() / "unr" -GPS_DIR.mkdir(exist_ok=True, parents=True) STATION_LLH_URL = "https://geodesy.unr.edu/NGLStationPages/llh.out" -STATION_LLH_FILE = str(GPS_DIR / "station_llh_all_{today}.csv") +STATION_LLH_FILENAME = "station_llh_all_{today}.csv" +STEPS_URL = "https://geodesy.unr.edu/NGLStationPages/steps.txt" +# Seconds before an HTTP download is abandoned (connect, read) +REQUEST_TIMEOUT = (10, 120) logger = logging.getLogger("geepers") @@ -48,7 +47,7 @@ def timeseries( frame: Literal["ENU", "XYZ"] = "ENU", start_date: str | None = None, end_date: str | None = None, - zero_by: Literal["mean", "start"] = "mean", + zero_by: Literal["mean", "start", "none"] = "mean", download_if_missing: bool = True, plate_fixed: bool = False, plate_name: str | None = None, @@ -86,7 +85,7 @@ def timeseries( msg = f"Unsupported frame: {frame}. Use 'ENU' or 'XYZ'" raise ValueError(msg) - station_id = station_id.upper() + station_id = validate_station_id(station_id.upper()) plate = None if plate_fixed and frame == "ENU": @@ -101,12 +100,12 @@ def timeseries( plate = plate_name else: plate = plates[0] - gps_data_file = GPS_DIR / plate / f"{station_id}.tenv3" + gps_data_file = self._cache_dir / plate / f"{station_id}.tenv3" else: if frame == "ENU": - gps_data_file = GPS_DIR / f"{station_id}.tenv3" + gps_data_file = self._cache_dir / f"{station_id}.tenv3" else: # frame in ("XYZ") - gps_data_file = GPS_DIR / f"{station_id}.txyz2" + gps_data_file = self._cache_dir / f"{station_id}.txyz2" if not gps_data_file.exists(): if download_if_missing: @@ -124,15 +123,7 @@ def timeseries( ) if frame == "ENU" and zero_by: - if zero_by.lower() == "mean": - mean_val = df[["east", "north", "up"]].mean() - df[["east", "north", "up"]] -= mean_val - elif zero_by.lower() == "start": - start_val = df[["east", "north", "up"]].iloc[:10].mean() - df[["east", "north", "up"]] -= start_val - else: - msg = "zero_by must be either 'mean' or 'start'" - raise ValueError(msg) + df = self._zero_data(df, zero_by, columns=["east", "north", "up"]) if frame == "ENU": StationObservationSchema.validate(df, lazy=True) @@ -149,8 +140,7 @@ def _read_station_data(self) -> gpd.GeoDataFrame: """ today = datetime.date.today().strftime("%Y%m%d") - filename = STATION_LLH_FILE.format(today=today) - lla_path = Path(filename) + lla_path = self._cache_dir / STATION_LLH_FILENAME.format(today=today) try: df = pd.read_csv(lla_path, sep=r"\s+", engine="c", header=None) @@ -194,27 +184,27 @@ def download_station_data( If None, uses the first plate from the UNR results. """ - station_id = station_id.upper() + station_id = validate_station_id(station_id.upper()) if frame == "ENU": if plate_fixed: if plate is None: plate = self._get_station_plates(station_id)[0] url = f"https://geodesy.unr.edu/gps_timeseries/tenv3/plates/{plate}/{station_id}.{plate}.tenv3" - filename = GPS_DIR / plate / f"{station_id}.tenv3" + filename = self._cache_dir / plate / f"{station_id}.tenv3" else: url = GPS_BASE_URL.format(station=station_id, reference=reference) # Hack to get around bad url structure url = url.replace("gps_timeseries/IGS14", "gps_timeseries") - filename = GPS_DIR / f"{station_id}.tenv3" + filename = self._cache_dir / f"{station_id}.tenv3" elif frame == "XYZ": url = f"https://geodesy.unr.edu/gps_timeseries/txyz/{reference}/{station_id}.txyz2" - filename = GPS_DIR / f"{station_id}.txyz2" + filename = self._cache_dir / f"{station_id}.txyz2" else: msg = "frame must be 'ENU' or 'XYZ'" raise ValueError(msg) - response = requests.get(url) + response = requests.get(url, timeout=REQUEST_TIMEOUT) response.raise_for_status() filename.parent.mkdir(parents=True, exist_ok=True) @@ -222,21 +212,29 @@ def download_station_data( logger.info(f"Saved {url} to {filename}") def _get_station_plates(self, station_id: str) -> list[str]: - """Get the tectonic plate for a given GPS station.""" - # A text file that gives the plate associated with each station is available: - url = "https://geodesy.unr.edu/gps_timeseries/Plates/sta_frames.txt" - # This directory also contains files for each frame ("plate_??.txt" where + """Get the tectonic plate(s) for a given GPS station.""" + plates = self._read_station_plates_table().get(station_id) + if plates is None: + msg = f"Failed to find {station_id} in the UNR station plates table" + raise ValueError(msg) + return plates + + @staticmethod + @cache + def _read_station_plates_table() -> dict[str, list[str]]: + """Download and parse the UNR station -> plates table (cached).""" + # A text file that gives the plate associated with each station. + # The directory also contains files for each frame ("plate_??.txt" where # ?? is the 2-character plate designation) that list the stations # associated with each plate. - response = requests.get(url) + url = "https://geodesy.unr.edu/gps_timeseries/Plates/sta_frames.txt" + response = requests.get(url, timeout=REQUEST_TIMEOUT) response.raise_for_status() + table: dict[str, list[str]] = {} for line in response.text.splitlines(): cur_id, *plates = line.split(" ") - if cur_id == station_id: - return plates - - msg = f"Failed to find {station_id} at {url}" - raise ValueError(msg) + table[cur_id] = plates + return table def _clean_gps_df( self, @@ -289,12 +287,78 @@ def _clean_gps_df( def _download_station_locations(self, filename: PathOrStr, url: str) -> None: """Download the station location file from the Nevada Geodetic Laboratory.""" - resp = requests.get(url) + resp = requests.get(url, timeout=REQUEST_TIMEOUT) resp.raise_for_status() with open(filename, "w") as f: f.write(resp.text) + def steps(self, station_ids: list[str] | None = None) -> pd.DataFrame: + """Fetch the UNR database of potential step epochs. + + Parses https://geodesy.unr.edu/NGLStationPages/steps.txt (format: + https://geodesy.unr.edu/NGLStationPages/steps_readme.txt), which + lists equipment changes (code 1) and earthquakes near the station + (code 2). + + Parameters + ---------- + station_ids : list of str, optional + Only return steps for these stations. + + Returns + ------- + pd.DataFrame + Columns: ``id``, ``date`` (parsed datetime), ``code``, + ``description``, plus for earthquake entries + ``threshold_distance``, ``distance_from_eq`` and + ``magnitude`` (NaN for equipment steps). + + """ + rows = self._read_steps_table() + df = rows.copy() + if station_ids is not None: + wanted = {s.upper() for s in station_ids} + df = df[df["id"].isin(wanted)].reset_index(drop=True) + return df + + @staticmethod + @cache + def _read_steps_table() -> pd.DataFrame: + """Download and parse the UNR steps file (cached).""" + response = requests.get(STEPS_URL, timeout=REQUEST_TIMEOUT) + response.raise_for_status() + + equipment, earthquakes = [], [] + for line in response.text.splitlines(): + parts = line.split() + if len(parts) < 4: + continue + if len(parts) > 5: # earthquake entries carry extra columns + earthquakes.append(parts[:7]) + else: + equipment.append(parts[:4]) + + df_eq = pd.DataFrame( + earthquakes, + columns=[ + "id", + "date", + "code", + "threshold_distance", + "distance_from_eq", + "magnitude", + "description", + ], + ) + df_env = pd.DataFrame(equipment, columns=["id", "date", "code", "description"]) + df = pd.concat([df_env, df_eq], ignore_index=True) + df["date"] = pd.to_datetime(df["date"], format="%y%b%d") + df["code"] = df["code"].astype(int) + for col in ("threshold_distance", "distance_from_eq", "magnitude"): + df[col] = pd.to_numeric(df[col], errors="coerce") + return df.sort_values(["id", "date"], ignore_index=True) + def get_global_station_list(self) -> list[str]: """Get the list of "processed" stations from UNR. diff --git a/src/geepers/gps_sources/unr_grid.py b/src/geepers/gps_sources/unr_grid.py index 788be0d..c45669a 100644 --- a/src/geepers/gps_sources/unr_grid.py +++ b/src/geepers/gps_sources/unr_grid.py @@ -2,10 +2,9 @@ from __future__ import annotations -from functools import lru_cache +from functools import cache, lru_cache from pathlib import Path -from sys import version -from typing import TYPE_CHECKING, Literal, Optional +from typing import TYPE_CHECKING, Literal import geopandas as gpd import pandas as pd @@ -13,40 +12,66 @@ from requests.adapters import HTTPAdapter, Retry from tqdm.contrib.concurrent import thread_map -from geepers.schemas import GridCellSchema, StationObservationSchema -from geepers.utils import decimal_year_to_datetime +from geepers.schemas import EPS, GridCellSchema, StationObservationSchema +from geepers.utils import decimal_years_to_datetimes from .base import BaseGpsSource +from .unr import REQUEST_TIMEOUT if TYPE_CHECKING: pass __all__ = ["UnrGridSource"] -VALID_VERSIONS = {"0.1", "0.2"} -DEFAULT_VERSION = "0.2" -LOOKUP_FILE_URL = "https://geodesy.unr.edu/grid_timeseries/Version{version}/grid_latlon_lookup.txt" +VALID_VERSIONS = {"0.1", "0.3"} +DEFAULT_VERSION: Literal["0.1", "0.3"] = "0.3" +# Grid points are only available for time-variable positions in version 0.1. +# Note: "time_contsant_gridded" is the actual (misspelled) directory name on +# the UNR server, not a typo introduced here. +GRIDDED_TYPE_DIRS = { + "constant": "time_contsant_gridded", + "variable": "time_variable_gridded", +} +LOOKUP_FILE_URL = ( + "https://geodesy.unr.edu/grid_timeseries/Version{version}/grid_latlon_lookup.txt" +) FILENAME_TEMPLATE = "{plate}/{grid_id:06d}_{plate}.tenv8" GRID_DATA_BASE_URL = ( - "https://geodesy.unr.edu/grid_timeseries/Version{version}/" - "time_variable_gridded/{filename}" + "https://geodesy.unr.edu/grid_timeseries/Version{version}/{gridded_dir}/{filename}" ) -# https://geodesy.unr.edu/grid_timeseries/time_variable_gridded/NA/000007_NA.tenv8 -# https://geodesy.unr.edu/grid_timeseries/time_variable_gridded/IGS14/000003_IGS14.tenv8 +# https://geodesy.unr.edu/grid_timeseries/Version0.3/time_variable_gridded/NA/000007_NA.tenv8 +# https://geodesy.unr.edu/grid_timeseries/Version0.3/time_contsant_gridded/IGS20/000003_IGS20.tenv8 class UnrGridSource(BaseGpsSource): """UNR Grid GPS data source for gridded time series data.""" - def __init__(self, - version: Literal["0.1", "0.2"] = "0.2", - cache_dir: Optional[str | Path] = None,): - """Initialize UNR Grid GPS data source.""" - + + def __init__( + self, + version: Literal["0.1", "0.3"] = DEFAULT_VERSION, + gridded_type: Literal["constant", "variable"] = "variable", + cache_dir: str | Path | None = None, + ): + """Initialize UNR Grid GPS data source. + + Parameters + ---------- + version : Literal["0.1", "0.3"], optional + Version of the UNR grid data to use. Default is "0.3". + gridded_type : Literal["constant", "variable"], optional + Whether to use the time-constant or time-variable gridded + products. Only available starting with version "0.3"; + version "0.1" only has time-variable data. Default is "variable". + cache_dir : str | Path, optional + Directory to cache downloaded data files. + + """ # Initialize BaseGpsSource super().__init__(cache_dir=cache_dir) # Store UNR version - self.version = version + self.version = version + self.gridded_type = gridded_type def timeseries( self, @@ -55,7 +80,7 @@ def timeseries( frame: Literal["ENU", "XYZ"] = "ENU", start_date: str | None = None, end_date: str | None = None, - zero_by: Literal["mean", "start"] = "mean", + zero_by: Literal["mean", "start", "none"] = "mean", download_if_missing: bool = True, *, plate: Literal["NA", "PA", "IGS14", "IGS20"] = "IGS14", @@ -96,12 +121,21 @@ def timeseries( # TODO: how to handle fetching/saving, vs using pandas to read... if download_if_missing: - local_file = self._download_file(station_id, plate=plate, - version=self.version) + local_file = self._download_file( + station_id, + plate=plate, + version=self.version, + gridded_type=self.gridded_type, + ) df = self.parse_data_file(local_file) else: - filename = FILENAME_TEMPLATE.format(plate=plate, grid_id=station_id) - uri = GRID_DATA_BASE_URL.format(version=version, filename=filename) + filename = FILENAME_TEMPLATE.format(plate=plate, grid_id=int(station_id)) + gridded_dir = GRIDDED_TYPE_DIRS[ + self.gridded_type if self.version != "0.1" else "variable" + ] + uri = GRID_DATA_BASE_URL.format( + version=self.version, gridded_dir=gridded_dir, filename=filename + ) df = self.parse_data_file(uri) df = self._filter_by_date(df, start_date, end_date) @@ -169,7 +203,8 @@ def _download_file( plate: Literal["NA", "PA", "IGS14", "IGS20"] = "IGS14", output_dir: Path | None = None, session: requests.Session | None = None, - version: Literal["0.1", "0.2"] = "0.2", + version: Literal["0.1", "0.3"] = DEFAULT_VERSION, + gridded_type: Literal["constant", "variable"] = "variable", ) -> Path: """Download ont .tenv8 data file. @@ -185,8 +220,12 @@ def _download_file( session : requests.Session A shared requests.Session object. Can be used for retrying. - version : Literal["0.1", "0.2"], optional + version : Literal["0.1", "0.3"], optional Version of the UNR grid data to download. + gridded_type : Literal["constant", "variable"], optional + Whether to download the time-constant or time-variable + gridded product. Version "0.1" only has time-variable data, + so this is ignored when version is "0.1". Returns ------- @@ -198,22 +237,25 @@ def _download_file( output_dir = self._cache_dir output_dir.mkdir(parents=True, exist_ok=True) - if plate == "IGS14" and version == "0.2": - # Warning: IGS14 plate is not available in UNR version 0.2. + if plate == "IGS14" and version == "0.3": + # Warning: IGS14 plate is not available in UNR version 0.3. # Using IGS20 instead. plate = "IGS20" - filename = FILENAME_TEMPLATE.format(plate=plate, grid_id=grid_id) + gridded_dir = GRIDDED_TYPE_DIRS[ + gridded_type if version != "0.1" else "variable" + ] + # Accept both int and zero-padded string ids ("000123" or 123) + filename = FILENAME_TEMPLATE.format(plate=plate, grid_id=int(grid_id)) url = GRID_DATA_BASE_URL.format( - version=version, - filename=filename + version=version, gridded_dir=gridded_dir, filename=filename ) dest = output_dir / url.rsplit("/", 1)[-1] if not dest.exists(): if session is None: - resp = requests.get(url) + resp = requests.get(url, timeout=REQUEST_TIMEOUT) else: - resp = session.get(url) + resp = session.get(url, timeout=REQUEST_TIMEOUT) resp.raise_for_status() with dest.open("wb") as f: f.write(resp.content) @@ -225,7 +267,8 @@ def download_data_files( plate: Literal["NA", "PA", "IGS14", "IGS20"] = "IGS14", max_workers: int = 8, output_dir: Path | None = None, - version: Literal["0.1", "0.2"] = "0.2", + version: Literal["0.1", "0.3"] = DEFAULT_VERSION, + gridded_type: Literal["constant", "variable"] = "variable", ) -> list[Path]: """Download .tenv8 data files in parallel, showing progress. @@ -241,8 +284,12 @@ def download_data_files( output_dir : Path | None, optional Directory to store downloaded data files. If None, the cache directory is used. - version : Literal["0.1", "0.2"], optional - Version of the UNR grid data to download. Default is "0.2". + version : Literal["0.1", "0.3"], optional + Version of the UNR grid data to download. Default is "0.3". + gridded_type : Literal["constant", "variable"], optional + Whether to download the time-constant or time-variable + gridded product. Ignored when version is "0.1", since that + version only has time-variable data. Returns ------- @@ -269,6 +316,7 @@ def download_data_files( session=s, desc="Downloading data files", version=version, + gridded_type=gridded_type, ) @staticmethod @@ -305,10 +353,11 @@ def parse_data_file(self, uri: str | Path) -> pd.DataFrame: DataFrame with columns validated against GPSUncertaintySchema. """ - df = self._read_data_file(uri) + # Copy so we never mutate the lru_cached DataFrame + df = self._read_data_file(uri).copy() - # Convert decimal year to datetime - df["date"] = df["decimal_year"].apply(decimal_year_to_datetime) + # Convert decimal year to datetime (vectorized) + df["date"] = decimal_years_to_datetimes(df["decimal_year"]) # Add placeholder correlation values (not in .tenv8 format) df["corr_en"] = 0.0 @@ -330,32 +379,37 @@ def parse_data_file(self, uri: str | Path) -> pd.DataFrame: "corr_nu", ] ] - # UNR Grid is in millimeters instead of meters: - df_out.loc[:, ["east", "north", "up"]] /= 1000 + # UNR Grid is in millimeters instead of meters (values and sigmas): + sigma_cols = ["sigma_east", "sigma_north", "sigma_up"] + df_out.loc[:, ["east", "north", "up", *sigma_cols]] /= 1000 + + # Time-constant gridded products can report exactly-zero sigmas, + # which the schema rejects; clamp to its minimum + df_out.loc[:, sigma_cols] = df_out[sigma_cols].clip(lower=EPS) return StationObservationSchema.validate(df_out, lazy=True) @staticmethod - @lru_cache(maxsize=None) - def _read_grid_file(version: str = "0.2") -> pd.DataFrame: + @cache + def _read_grid_file(version: str = DEFAULT_VERSION) -> pd.DataFrame: """Download and cache the UNR grid latitude/longitude lookup table.""" if version not in VALID_VERSIONS: msg = ( f"Unsupported version '{version}'. " f"Valid options are: {', '.join(VALID_VERSIONS)}." - ) - raise ValueError(msg) - + ) + raise ValueError(msg) + url = LOOKUP_FILE_URL.format(version=version) df = pd.read_csv( url, sep=r"\s+", names=["grid_point", "longitude", "latitude"], ) - if version == "0.2": - # Version 0.2 maps latitudes from 0 to 360 + if version == "0.3": + # Version 0.3 maps longitudes from 0 to 360 # convert to -180 to 180, to be consistent with version 0.1 - df['longitude'] = ((df['longitude'] + 180) % 360) - 180 + df["longitude"] = ((df["longitude"] + 180) % 360) - 180 return df.set_index("grid_point") diff --git a/src/geepers/utils.py b/src/geepers/utils.py index 30d7caf..fd9dd4f 100644 --- a/src/geepers/utils.py +++ b/src/geepers/utils.py @@ -7,6 +7,7 @@ import geopandas as gpd import numpy as np +import pandas as pd from ._types import DateOrDatetime @@ -78,7 +79,7 @@ def get_cache_dir(force_posix=False, app_name="geepers") -> Path: return Path("~/Library/Application Support").expanduser() / app_name else: base_path = Path( - os.environ.get("XDG_CONFIG_HOME", Path("~/.cache").expanduser()) + os.environ.get("XDG_CACHE_HOME", Path("~/.cache").expanduser()) ) return base_path / app_name @@ -201,3 +202,23 @@ def decimal_year_to_datetime(decimal_year: float) -> datetime.datetime: return datetime.datetime(1990, 1, 1) + datetime.timedelta( seconds=(decimal_year - start_year) * seconds_per_year ) + + +def decimal_years_to_datetimes(decimal_years) -> pd.DatetimeIndex: + """Vectorized version of `decimal_year_to_datetime`. + + Parameters + ---------- + decimal_years : array-like of float + Years expressed as decimals (e.g., [2014.5, 2014.5027]). + + Returns + ------- + pd.DatetimeIndex + Corresponding calendar datetimes, using the same 365.25-day-year + convention as `decimal_year_to_datetime`. + + """ + dy = np.asarray(decimal_years, dtype=float) + seconds = (dy - 1990.0) * 365.25 * 24 * 3600 + return pd.Timestamp("1990-01-01") + pd.to_timedelta(seconds, unit="s") From abba1eec51826a321ca67d17e8e15aa0cd594e0a Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Mon, 13 Jul 2026 12:18:11 -0700 Subject: [PATCH 2/4] Make test_main row assertion order-independent The spatial-filter method (cx bbox-slice) changes station block order in combined_data.csv, so HLNA's row is no longer at index 0. Compare row content with check_names=False instead of asserting the absolute label. Co-Authored-By: Claude Opus 4.8 --- tests/test_workflows.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/tests/test_workflows.py b/tests/test_workflows.py index 8cc7fde..f18b313 100644 --- a/tests/test_workflows.py +++ b/tests/test_workflows.py @@ -43,8 +43,12 @@ def test_main(tmp_path, monkeypatch): "measurement": "los_gps", "value": 0.0010796103397325, } + # Compare row content independent of its absolute position in the CSV + # (station block order depends on the spatial-filter method). pd.testing.assert_series_equal( - df[df.id == "HLNA"].iloc[0], pd.Series(expected_entry, name=0) + df[df.id == "HLNA"].iloc[0], + pd.Series(expected_entry, name=0), + check_names=False, ) From b30906aa4d173acd9343e993f307d9c76c9b1dfd Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Mon, 13 Jul 2026 15:18:31 -0700 Subject: [PATCH 3/4] Deploy mkdocs docs at /docs/ alongside the viewer deploy-pages.sh now builds the mkdocs site into docs/ (viewer stays at the Pages root); .nojekyll at root covers the docs assets. Viewer links to Docs. Co-Authored-By: Claude Opus 4.8 --- scripts/browse_unr_grid.html | 1 + scripts/deploy-pages.sh | 16 +++++++++++++++- 2 files changed, 16 insertions(+), 1 deletion(-) diff --git a/scripts/browse_unr_grid.html b/scripts/browse_unr_grid.html index 7e4ea31..ef78d13 100644 --- a/scripts/browse_unr_grid.html +++ b/scripts/browse_unr_grid.html @@ -321,6 +321,7 @@

OPERA UNR Grid Viewer Gridded GPS time series by the Nevada Geodetic Laboratory (UNR), funded by the JPL-led OPERA project. Viewer: JPL. + · Docs

diff --git a/scripts/deploy-pages.sh b/scripts/deploy-pages.sh index 03d29a8..42d672f 100755 --- a/scripts/deploy-pages.sh +++ b/scripts/deploy-pages.sh @@ -48,7 +48,20 @@ PY # 2. Static assets. cp "$DATA_SRC" "$WORK/$DATA_DEST" [ -f "$BOUNDARIES" ] && cp "$BOUNDARIES" "$WORK/PB2002_boundaries.json" +# .nojekyll (at root) disables Jekyll for the whole site, so the mkdocs +# assets under docs/ (e.g. _mkdocstrings.css) are served too. : > "$WORK/.nojekyll" + +# 2b. mkdocs docs at /docs/ (viewer stays at the site root). Built only when +# mkdocs is available (run from an env with the docs deps + geepers importable); +# skipped with a warning otherwise so a viewer-only deploy still works. +REPO_ROOT="$(cd "$SCRIPTS_DIR/.." && pwd)" +if command -v mkdocs >/dev/null 2>&1 && [ -f "$REPO_ROOT/mkdocs.yml" ]; then + echo "Building docs → docs/ …" + ( cd "$REPO_ROOT" && PYTHONPATH=src mkdocs build --quiet --site-dir "$WORK/docs" ) +else + echo "warning: mkdocs not found; deploying viewer only (no docs/)." >&2 +fi cat > "$WORK/README.md" <\` (host must allow CORS + ranges). EOF @@ -67,7 +81,7 @@ EOF cd "$WORK" git init -q -b gh-pages git add -A -git commit -q -m "Deploy OPERA UNR grid viewer ($(basename "$DATA_SRC"), ${size_mb} MB)" +git commit -q -m "Deploy OPERA UNR grid viewer ($(basename "$DATA_SRC"), ${size_mb} MB) + docs" git remote add origin "$REMOTE" git push -f origin gh-pages From f3a3ad509350379a677176fdc8832896a827d368 Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Mon, 13 Jul 2026 15:19:58 -0700 Subject: [PATCH 4/4] Inject viewer Docs link at deploy time only (no dead link when run locally) Co-Authored-By: Claude Opus 4.8 --- scripts/browse_unr_grid.html | 1 - scripts/deploy-pages.sh | 15 +++++++++++++-- 2 files changed, 13 insertions(+), 3 deletions(-) diff --git a/scripts/browse_unr_grid.html b/scripts/browse_unr_grid.html index ef78d13..7e4ea31 100644 --- a/scripts/browse_unr_grid.html +++ b/scripts/browse_unr_grid.html @@ -321,7 +321,6 @@

OPERA UNR Grid Viewer Gridded GPS time series by the Nevada Geodetic Laboratory (UNR), funded by the JPL-led OPERA project. Viewer: JPL. - · Docs

diff --git a/scripts/deploy-pages.sh b/scripts/deploy-pages.sh index 42d672f..311970c 100755 --- a/scripts/deploy-pages.sh +++ b/scripts/deploy-pages.sh @@ -34,7 +34,9 @@ fi WORK="$(mktemp -d)" trap 'rm -rf "$WORK"' EXIT -# 1. Viewer: point the default DATA_URL at the hosted (renamed) dataset. +# 1. Viewer: point the default DATA_URL at the hosted (renamed) dataset, and +# inject a "Docs" link (deploy-only, so running the HTML locally has no +# dead docs/ link). python3 - "$SCRIPTS_DIR/browse_unr_grid.html" "$WORK/index.html" "$DATA_DEST" <<'PY' import sys src, dst, data = sys.argv[1:4] @@ -42,7 +44,16 @@ html = open(src).read() old = "const DATA_URL = params.get('data') || 'unr_grid.parquet';" new = f"const DATA_URL = params.get('data') || '{data}';" assert html.count(old) == 1, "could not find DATA_URL default in viewer" -open(dst, "w").write(html.replace(old, new)) +html = html.replace(old, new) +# Deploy-only Docs link (present only when docs/ is published alongside). +credit = "funded by the JPL-led OPERA project. Viewer: JPL." +if credit in html: + html = html.replace( + credit, + credit + '\n · ' + 'Docs', + ) +open(dst, "w").write(html) PY # 2. Static assets.