diff --git a/SPEC.md b/SPEC.md index b808d66..1343f87 100644 --- a/SPEC.md +++ b/SPEC.md @@ -94,11 +94,11 @@ Cloud Scheduler (cron) GCS bucket (public read) │ - GDAL /vsigs/ (GeoJSON) · PyArrow gs:// (GeoParquet) + copied into the image at build (Cloud Build) │ Cloud Run Service ← pygeoapi (always-on) - (OGR provider: GeoJSON products; - Parquet provider: timeseries products) + (Parquet provider: point products; + OGR provider: polygon products, GeoJSON) │ HTTP clients (OGC API - Features) @@ -106,11 +106,12 @@ Cloud Scheduler (cron) **Two Cloud Run deployments:** - **Cloud Run Job** (Dagster, stateless): triggered by Cloud Scheduler, runs pipeline per product -- **Cloud Run Service** (pygeoapi, always-on): serves OGC API - Features from GCS +- **Cloud Run Service** (pygeoapi, always-on): serves OGC API - Features from product + files baked into its image -No PostgreSQL required. GCS is authoritative store. pygeoapi reads product files directly -from GCS — GeoJSON via GDAL's `/vsigs/` virtual filesystem, timeseries GeoParquet via -PyArrow's `gs://` filesystem — no proxy, no DB, no sync step (§6). +No PostgreSQL required. GCS is authoritative store. The pygeoapi image build copies each +product's current file from GCS into the image, so the service reads local disk — no +proxy, no DB — and new data is served after the next image build (§6.3, §6.6). > GeoServer/PostGIS persister in `backend/persisters/geoserver.py` is unchanged and > still usable via CLI `--output-format geoserver`. Not part of this pipeline. @@ -186,7 +187,8 @@ temporal filtering natively without custom code. `backend/persisters/ogc_features.py` → `OGCFeaturesPersister` - `dump_summary_collection(path, records, meta)` — §4.2 format - `dump_timeseries_collection(path, site_records, timeseries_records, meta)` — §4.3 format -- Writes local `.geojson` file; Dagster GCS resource handles upload +- Writes local `.geojson` file; Dagster GCS resource handles upload, and also writes a + GeoParquet copy (`latest.parquet`, see §6.1) via `geodataframe.collection_to_geoparquet` --- @@ -207,6 +209,7 @@ products_config ← loads products.yaml at startup ▼ gcs_upload ← gs://die-products/{product_id}/{YYYY-MM-DD}.geojson gs://die-products/{product_id}/latest.geojson (overwrite) + gs://die-products/{product_id}/latest.parquet (overwrite) ``` ### §5.2 Configurable Products (`orchestration/config/products.yaml`) @@ -264,32 +267,56 @@ Structure supports later migration to persistent Dagster daemon (change Cloud Ru ### §6.1 Role pygeoapi serves the GCS-stored product files as OGC API - Features collections. No DB. -The provider depends on the product's `output_type`: +Each product is served from one of two files, chosen and copied into the image when it +is built (§6.3): -| Products | Provider | Source | +| When | Provider | Source (in the image; copied from `gs://{bucket}/products/{id}/`) | |---|---|---| -| `ogc_timeseries` | **Parquet** (`time_field: datetime`) | `gs://{bucket}/products/{id}/latest.parquet` | -| all others | **OGR** (`source_type: GeoJSON`) | `/vsigs/{bucket}/products/{id}/latest.geojson` | - -Timeseries products use Parquet because the OGR provider ignores `datetime` filters and -parses the whole GeoJSON on every request (~9–10 s per request for the 733 MB -`nm_waterlevels_timeseries` file in local testing, vs. under 1 s from GeoParquet). - -The timeseries GeoParquet file must have: -- GeoParquet `geo` metadata with WKB geometry and a **bbox covering column** (required for - `bbox` queries; geopandas `to_parquet(..., write_covering_bbox=True)`) -- `datetime` as a **tz-aware UTC** timestamp -- `string` rather than `large_string` columns (the provider can't describe `large_string` - fields) - -The pipeline does not write `latest.parquet` yet; it is uploaded manually. It must exist -in the bucket before pygeoapi starts — the entrypoint fails if any collection can't open. +| `latest.parquet` is usable | **Parquet** (`id_field: feature_id`; `time_field: datetime` for `ogc_timeseries`) | `/data/products/{id}/latest.parquet` | +| otherwise | **OGR** (`source_type: GeoJSON`) | `/data/products/{id}/latest.geojson` | + +Parquet is preferred because the OGR provider parses the whole GeoJSON on every request +(~9–10 s for the 733 MB `nm_waterlevels_timeseries` file, ~0.1 s per MB from GCS for the +rest, vs. well under 1 s from GeoParquet) and ignores `datetime` filters. + +**Writing `latest.parquet`.** `GCSResource.upload_product` writes it after every GeoJSON +upload, stamped with the GeoJSON's content hash and the converter's +`PARQUET_FORMAT_VERSION`. When the GeoJSON is unchanged it still writes the Parquet if it +is missing or either stamp differs (backfill; bump the version when the converter's output +changes). A conversion failure +does not fail the product: the Parquet is deleted so a stale copy can't be served, and the +asset's `parquet_status` / `parquet_error` metadata record it. The file has: +- GeoParquet 1.1 `geo` metadata with WKB geometry and a **bbox covering column** +- `feature_id`: the feature's top-level GeoJSON id, made unique (repeats get `:2`, `:3`, …; + timeseries ids are `source:site:date`, which repeats for several readings a day) +- `datetime` (when present) as a **tz-aware UTC** `timestamp[ms]`, rows sorted by `id`, + `source`, then `datetime`, so a one-well query (`?id=`) reads one or two row groups +- `string` rather than `large_string` columns; dict/list values JSON-encoded; columns + mixing bool/number/string values stored as strings + +**Usable** means: GeoParquet metadata with a bbox covering column, **point geometry only**, +an id column (`feature_id`, else `id`), and for `ogc_timeseries` a tz-aware `datetime`. +Polygon layers stay on GeoJSON because the Parquet provider answers `bbox` by containment +(the feature's bbox fully inside the query box), not intersection, and would drop +polygons crossing the box edge. + +The image patches the Parquet provider (`orchestration/pygeoapi/patch_pygeoapi.py`, run by +the Dockerfile against the pinned base image; the build fails if the patched source +changes): +- **Paging:** a filtered page returned up to `limit + 1` rows plus earlier batches, so + pages overlapped. Every page but the last now has exactly `limit` rows. +- **Date-only `datetime`** bounds (`2020-01-01/2020-12-31`) returned 500 against the + tz-aware column. They're now read as UTC. +- **No `bbox`:** responses carried a collection `bbox`, which was `[NaN, NaN, NaN, NaN]` + (invalid JSON) for empty results, and a `bbox` on every feature. Both are now left + out, matching the OGR collections. +- **True `numberMatched`:** items pages report the total count of matching rows (as + `resulttype=hits` does), not `offset` + rows returned + 1. That costs ~20–110 ms of + server time per request. Known provider limitations: -- Parquet: date-only `datetime` queries (`2020-01-01/2020-12-31`) return 500 against a - tz-aware column — full RFC 3339 values (`2020-01-01T00:00:00Z/..`) work. - `numberMatched` is wrong. -- OGR: no `numberReturned`; DateTime properties are returned as `YYYY/MM/DD HH:MM:SS+00`. +- OGR: no `numberMatched` or `numberReturned`; DateTime properties are returned as + `YYYY/MM/DD HH:MM:SS+00`. ``` GET /collections @@ -312,6 +339,7 @@ server: language: en-US cors: true pretty_print: false + gzip: true # Cloud Run doesn't compress limits: default_items: 500 max_items: 10000 # NewWeaver requests limit=10000 @@ -324,8 +352,12 @@ logging: metadata: identification: - title: NM Unified Water Data - description: OGC API - Features for New Mexico water data + title: New Mexico Unified Water Data + description: >- + This tool integrates groundwater-level and water-quality data for New Mexico + from state and federal sources, including NMBGMR, USGS, and the Water Quality Portal. + Its collections provide per-well summaries, trends, time series, and derived + water-quality indicators through the OGC API - Features standard. keywords: - water - groundwater @@ -333,7 +365,7 @@ metadata: - NMBGMR keywords_type: theme terms_of_service: https://creativecommons.org/licenses/by/4.0/ - url: https://waterdata.nmt.edu + url: ${PYGEOAPI_SERVER_URL} # the live Cloud Run URL, set in cloudbuild.yaml license: name: CC-BY 4.0 url: https://creativecommons.org/licenses/by/4.0/ @@ -347,6 +379,7 @@ metadata: stateorprovince: New Mexico postalcode: "87801" country: USA + email: data-services-nmbg@nmt.edu resources: {% for product in products %} @@ -372,22 +405,22 @@ resources: end: null {% endif %} providers: -{% if product.output_type == 'ogc_timeseries' %} - # GeoParquet via the Parquet provider: the OGR provider ignores datetime - # filters and parses the whole GeoJSON on every request. The file must be - # GeoParquet with a bbox covering column and a tz-aware UTC datetime. +{% if product.id in parquet %} + # Chosen by generate_config.py when latest.parquet is usable. - type: feature name: Parquet data: - source: gs://{{ gcs_bucket }}/products/{{ product.id }}/latest.parquet - id_field: id + source: {{ parquet_root }}/{{ product.id }}/latest.parquet + id_field: {{ parquet[product.id] }} +{% if product.output_type == 'ogc_timeseries' %} time_field: datetime +{% endif %} {% else %} - type: feature name: OGR data: source_type: GeoJSON - source: /vsigs/{{ gcs_bucket }}/products/{{ product.id }}/latest.geojson + source: {{ geojson_root }}/{{ product.id }}/latest.geojson source_options: GDAL_HTTP_UNSAFESSL: "NO" gdal_ogr_options: @@ -409,50 +442,42 @@ basename). ### §6.3 Config Generation (`orchestration/pygeoapi/generate_config.py`) -```python -def generate(products_path: Path, template_path: Path, output_path: Path) -> None: - products_config = yaml.safe_load(products_path.read_text()) - env = Environment( - loader=FileSystemLoader(str(template_path.parent)), - keep_trailing_newline=True, - ) - tmpl = env.get_template(template_path.name) - rendered = tmpl.render( - products=products_config["products"], - gcs_bucket=products_config["gcs_bucket"], - ) - - # Sanity check: every product must produce a GCS source entry - for product in products_config["products"]: - pid = product["id"] - bucket = products_config["gcs_bucket"] - if product.get("output_type") == "ogc_timeseries": - expected = f"gs://{bucket}/products/{pid}/latest.parquet" - else: - expected = f"/vsigs/{bucket}/products/{pid}/latest.geojson" - assert expected in rendered, ( - f"§V violated: provider for '{pid}' missing {expected} in generated config" - ) - - output_path.parent.mkdir(parents=True, exist_ok=True) - output_path.write_text(rendered) - print(f"Generated {output_path} ({len(products_config['products'])} collections)") -``` +Renders `config.yml.j2` from `products.yaml` at **image build time** — baked into the +image, not runtime. It reads each product's `latest.parquet` footer from GCS and picks +the Parquet provider when the file is usable (§6.1), otherwise OGR + GeoJSON; the build +log lists each product's choice and reason. A missing or unreadable file falls back to +GeoJSON; any other error (permissions, network) fails the build rather than silently +serving everything as GeoJSON. + +- `--bake-to /data/products` (Cloud Build, `BAKE_PRODUCTS=1`): also copies the chosen + file into the image and points the config at the copy. A product with neither a usable + `latest.parquet` nor a `latest.geojson` fails the build, so the previous revision keeps + serving. New data is served after the next build (§6.6). +- `--check-parquet` (`CHECK_PARQUET=1`): same choice, but the config reads GCS directly + (OGR via `/vsigs/`, Parquet via `gs://`). +- Neither (e.g. a plain local `docker build`): every product is GeoJSON read from GCS. -Run at Docker build time — baked into the image, not runtime. +`--source-root DIR` replaces GCS with a local directory — local testing only. ### §6.4 GCS Auth +With baked products (the Cloud Build default) the service reads no GCS at runtime; only +the build does. The rest of this section applies to images built without `BAKE_PRODUCTS`. + On Cloud Run, both GDAL `/vsigs/` and PyArrow's `gs://` filesystem use the service account's ADC. No credentials file needed. GDAL only uses the metadata server when it detects GCE, which it can't on Cloud Run — `CPL_MACHINE_IS_GCE=YES` must be set or every -OGR collection fails with `No valid GCS credentials found`. PyArrow needs no setting. Require the pygeoapi Cloud Run -Service account to have `roles/storage.objectViewer` on the products bucket -(`gcs_bucket` in `products.yaml`). +OGR collection fails with `No valid GCS credentials found`. PyArrow needs no setting. +Require the pygeoapi Cloud Run service account to have `roles/storage.objectViewer` on the products bucket +(`gcs_bucket` in `products.yaml`). The Cloud Build account needs the same role for the +build-time check and copy (§6.3). -For local dev, render the config and rewrite GCS paths to a mounted folder instead: +For local dev, render the config inside the image against a mounted products folder: ```bash -sed -i '' -E 's#(/vsigs/|gs://)dataservices-die-products#/data#' local.config.yml +docker run --rm --entrypoint /venv/bin/python3 \ + -v ~/pygeoapi-test/products:/data/products:ro -v ~/pygeoapi-test:/out die-pygeoapi \ + /tmp/generate_config.py --products /tmp/products.yaml --template /tmp/config.yml.j2 \ + --output /out/local.config.yml --check-parquet --source-root /data/products ``` ### §6.5 Dockerfile (`orchestration/pygeoapi/Dockerfile`) @@ -460,10 +485,11 @@ sed -i '' -E 's#(/vsigs/|gs://)dataservices-die-products#/data#' local.config.ym Build context is `orchestration/` so `config/products.yaml` can be copied. ```dockerfile -FROM geopython/pygeoapi:latest +# Pinned for patch_pygeoapi.py: pygeoapi 0.25.dev0 (:latest on 2026-09-24). +FROM geopython/pygeoapi@sha256:f3dd50a56f870d80df67416b6b07bf4d92aefa78b7e6f697564288f1cdf1df3f -# pygeoapi base image includes GDAL with /vsigs/ GCS support. -# Auth uses Application Default Credentials — no key file needed on Cloud Run. +# Unbaked images read GCS directly (GDAL /vsigs/, Application Default +# Credentials). WORKDIR /pygeoapi @@ -473,6 +499,10 @@ WORKDIR /pygeoapi RUN /venv/bin/python3 -m pip install --no-cache-dir --quiet \ pyarrow==25.0.1 geopandas==1.1.4 +# Fix Parquet provider bugs; fails the build if the patched source changed. +COPY pygeoapi/patch_pygeoapi.py /tmp/patch_pygeoapi.py +RUN /venv/bin/python3 /tmp/patch_pygeoapi.py + # Copy generation inputs # Build context is orchestration/ (see cloudbuild.yaml) COPY pygeoapi/config.yml.j2 /tmp/config.yml.j2 @@ -482,10 +512,19 @@ COPY config/products.yaml /tmp/products.yaml # Bake config into image at build time. # §V: config generated from products.yaml, not hand-edited. # jinja2 and PyYAML are pygeoapi dependencies, already in /venv. +# BAKE_PRODUCTS=1 (cloudbuild.yaml): copy each product's file from GCS into +# /data/products and serve it locally; needs GCS access at build time. Use +# --no-cache locally to re-download. +# CHECK_PARQUET=1: choose Parquet vs GeoJSON the same way, but read GCS. +# Neither: GeoJSON from GCS. +ARG BAKE_PRODUCTS=0 +ARG CHECK_PARQUET=0 RUN /venv/bin/python3 /tmp/generate_config.py \ --products /tmp/products.yaml \ --template /tmp/config.yml.j2 \ - --output /pygeoapi/local.config.yml + --output /pygeoapi/local.config.yml \ + $([ "$CHECK_PARQUET" = "1" ] && echo --check-parquet) \ + $([ "$BAKE_PRODUCTS" = "1" ] && echo --bake-to /data/products) EXPOSE 80 @@ -493,14 +532,33 @@ EXPOSE 80 ENV PYGEOAPI_CONFIG=/pygeoapi/local.config.yml ``` +Cloud Build (`orchestration/pygeoapi/cloudbuild.yaml`) builds with `--network=cloudbuild` and +`--build-arg=BAKE_PRODUCTS=1` so the build can check and copy the products from GCS (§6.3). +Each Cloud Build runs a fresh Docker daemon, so no layer cache serves stale data; locally, +build with `--no-cache` to re-download. + Cloud Run Service (`orchestration/pygeoapi/cloudbuild.yaml`): - `PYGEOAPI_SERVER_URL` — `https://die-pygeoapi-$PROJECT_NUMBER.$_REGION.run.app` -- `CPL_MACHINE_IS_GCE=YES` — lets GDAL `/vsigs/` use Cloud Run credentials (§6.4) +- `CPL_MACHINE_IS_GCE=YES` — lets GDAL `/vsigs/` use Cloud Run credentials (§6.4); unused + with baked products, kept for images built without them - Port: 80 - Service account: `die-pygeoapi@$PROJECT_ID.iam.gserviceaccount.com` — needs `roles/storage.objectViewer` on the products bucket; Cloud Build's account needs `roles/iam.serviceAccountUser` on it -- Memory: 2Gi — the Parquet provider peaked at ~1.8 GB across 4 gunicorn workers locally +- `WSGI_WORKERS=2` — gunicorn workers (base image default 4) +- Memory: 2Gi — with 4 workers a sustained load of 10,000-row pages pushed the local + container to its 2 GiB cap; 2 workers peaked at ~1.8 GiB with the same timings + +### §6.6 Rebuild After a Pipeline Run + +Baked products (§6.3) are only as fresh as the image: a pipeline run that publishes new +data (a combine asset with `skipped_unchanged: false` or `parquet_status: written`) is +not served until pygeoapi is rebuilt and redeployed with +`orchestration/pygeoapi/cloudbuild.yaml`. **For now this rebuild is manual**: run the +build after the scheduled product runs finish, and not while any are still running, so +the image doesn't copy a half-updated set of products. Automating it (a Cloud Build +trigger started from Dagster after product runs publish new data) is planned; see +`docs/pygeoapi-performance-plan.md` T5. --- @@ -1049,8 +1107,8 @@ feature/composition-refactor ← branch off main after §T.9 merged - `die` CLI behavior unchanged after uv migration - All existing tests pass under `uv run pytest` - pygeoapi config MUST be generated from `products.yaml` — never hand-edited -- pygeoapi OGR provider MUST use `/vsigs/` path (GCS), never local filesystem path -- pygeoapi Parquet provider (`ogc_timeseries` products) MUST use `gs://` path (GCS), never local filesystem path +- pygeoapi product files MUST come from GCS (`gs://{bucket}/products/{id}/latest.*`): read + directly, or copied into the image at build time (§6.3) — never committed or hand-copied - No database introduced in orchestration pipeline — GCS is sole store - `latest.geojson` MUST be overwritten atomically (upload to tmp key, then copy/rename) diff --git a/backend/persisters/geodataframe.py b/backend/persisters/geodataframe.py index c39f523..0366ca5 100644 --- a/backend/persisters/geodataframe.py +++ b/backend/persisters/geodataframe.py @@ -24,6 +24,7 @@ import io import json +from datetime import datetime, timezone from pathlib import Path import geopandas as gpd @@ -292,3 +293,112 @@ def geojson_to_geopackage(geojson_path, layer_name: str, out_dir) -> tuple: gpkg_path = Path(out_dir) / f"{layer_name}.gpkg" bbox = write_geopackage(gdf, str(gpkg_path), layer_name) return gpkg_path, bbox + + +# --------------------------------------------------------------------------- +# GeoParquet output for pygeoapi's Parquet provider (SPEC §6.1). +# --------------------------------------------------------------------------- + +PARQUET_ID_FIELD = "feature_id" +PARQUET_TIME_FIELD = "datetime" +# Sort keys before datetime, so an `id=` filter skips row groups. id comes +# first because clients filter on it alone. +PARQUET_SITE_FIELDS = ("id", "source") +# Bump when the output changes, so existing latest.parquet files get rewritten. +PARQUET_FORMAT_VERSION = "2" + + +def _utc_timestamp_ms(value): + """ISO-8601 string -> UTC datetime (naive values taken as UTC), or None. + Parsed in Python: pandas ns timestamps overflow past 2262, and Arrow + rejects offset-less strings for a tz-aware type.""" + if _is_null(value): + return None + try: + dt = datetime.fromisoformat(str(value)) + except ValueError: + return None + return dt.replace(tzinfo=timezone.utc) if dt.tzinfo is None else dt.astimezone(timezone.utc) + + +def _unique_feature_ids(features: list[dict]) -> list[str]: + """Feature ids as strings, with repeats suffixed ``:2``, ``:3``, … so + ``/items/{id}`` finds one feature. A missing id becomes the position.""" + seen: dict[str, int] = {} + out = [] + for i, f in enumerate(features, start=1): + fid = str(f["id"]) if f.get("id") is not None else str(i) + n = seen.get(fid, 0) + 1 + seen[fid] = n + out.append(fid if n == 1 else f"{fid}:{n}") + return out + + +def _scalar_kind(v) -> str: + if isinstance(v, bool): + return "bool" + if isinstance(v, (int, float)): + return "number" + return type(v).__name__ + + +def _stringify_multi_type_column(gdf: pd.DataFrame, col: str) -> None: + """Cast *col* to str in place if it mixes bool/number/str values + (e.g. floats with 'ND'); Arrow needs one type per column.""" + kinds = {_scalar_kind(v) for v in gdf[col] if not _is_null(v)} + if len(kinds) > 1: + gdf[col] = [None if _is_null(v) else str(v) for v in gdf[col]] + + +def collection_to_geoparquet(collection: dict, out_path) -> int: + """Write a parsed product GeoJSON *collection* to GeoParquet 1.1 (with a + bbox covering column) at *out_path*. Returns the row count. + + Adds a unique ``feature_id``; stores ``datetime`` as UTC ``timestamp[ms]`` + and sorts by site then time; JSON-encodes dict/list values; stringifies + mixed-type columns.""" + import pyarrow as pa + import pyarrow.parquet as pq + + features = collection.get("features", []) + gdf = gpd.GeoDataFrame.from_features(features, crs="EPSG:4326") + gdf[PARQUET_ID_FIELD] = _unique_feature_ids(features) + + for col in gdf.columns: + if col == "geometry": + continue + if any(isinstance(v, (dict, list)) for v in gdf[col]): + gdf[col] = [ + json.dumps(v) if isinstance(v, (dict, list)) else v for v in gdf[col] + ] + _stringify_multi_type_column(gdf, col) + + if PARQUET_TIME_FIELD in gdf.columns: + ts = pa.array( + [_utc_timestamp_ms(v) for v in gdf[PARQUET_TIME_FIELD]], + type=pa.timestamp("ms", tz="UTC"), + ) + gdf[PARQUET_TIME_FIELD] = pd.arrays.ArrowExtensionArray(ts) + keys = [c for c in PARQUET_SITE_FIELDS if c in gdf.columns] + gdf = gdf.sort_values(keys + [PARQUET_TIME_FIELD], kind="stable") + + gdf.to_parquet( + out_path, + index=False, + schema_version="1.1.0", + write_covering_bbox=True, + row_group_size=100_000, + ) + + # pandas 3 writes large_string, which the provider can't read; downcast. + table = pq.read_table(out_path) + if any(pa.types.is_large_string(f.type) for f in table.schema): + schema = pa.schema( + [ + f.with_type(pa.string()) if pa.types.is_large_string(f.type) else f + for f in table.schema + ], + metadata=table.schema.metadata, + ) + pq.write_table(table.cast(schema), out_path, row_group_size=100_000) + return table.num_rows diff --git a/docs/pygeoapi-performance-plan.md b/docs/pygeoapi-performance-plan.md new file mode 100644 index 0000000..36fec66 --- /dev/null +++ b/docs/pygeoapi-performance-plan.md @@ -0,0 +1,302 @@ +# Plan: pygeoapi performance and correctness (post-standup) + +Status: T1–T4 done on DIE-363 and deployed (2026-09-24); T5 deferred (rebuild is manual +for now); T6 partly (format version); T7 partly (per-feature bbox). See "Progress" below. +Picks up from the pygeoapi standup (PR #137, branch +`DIE-359-deploy-pygeo-api-on-cloud-run`) and the per-product GeoParquet work (branch +`DIE-363-update-the-layer-generation-code-to-land-in-pygeo-api-instead-of-geo-server`, +4 commits on top of #137, not yet pushed). Read `SPEC.md` §6 first; it describes the +current design. + +## Progress (2026-09-24) + +Measured locally: patched image, files on local disk, 4 workers, 2 GiB. + +| Task | State | Result | +|---|---|---| +| T1 gzip | done | `id=9033` page: 3.85 MB → 83 KB on the wire; unfiltered `limit=5000`: 4.25 MB → 303 KB | +| T2 pin + patch | done | pages are exactly 5,000 rows, no overlap, 48,809 distinct ids; date-only `datetime` returns 200 | +| T3 sort | done, changed | sorted by **`id`, `source`**, then `datetime` (not `source` first: with source first the `id` ranges of row groups overlap and well 9033 hit 4 of 9 groups; id first, 1–2). Whole series 1.3–1.4 s vs 1.6–2.5 s time-sorted | +| T4 bake | done | `--bake-to`; 31/31 collections pass items/bbox/by-id from baked files; 88 MB baked; `limit=10` in 13–35 ms | +| T5 rebuild | deferred | Not automated yet. Rebuild manually after product runs publish new data (SPEC §6.6) | +| T6 | partial | `PARQUET_FORMAT_VERSION` done (now `"2"`, so every existing `latest.parquet` is rewritten on the next run). Others open | +| T7 | partial | per-feature `bbox` removed (with the empty-result fix below). The `bbox` covering column is still a property | + +**Live results** (baked build with the locally built Parquet files uploaded by hand, 2 +workers): the whole well-9033 series in 5,000-row pages takes 9.6–9.7 s (GeoServer ~10 s; +before, ~26 s), with 48,809 distinct rows and no overlap. With 10,000-row pages it takes +7.9–8.0 s: each request has ~0.4 s of fixed cost, so fewer, larger pages are faster. +Most `limit=10` requests take 0.17–0.22 s total, of which ~0.13 s is network. Cloud Run +memory stayed at 26–29% of 2 GiB, and the latency breakdown is almost all user execution. + +**Client feedback, fixed after the first live test** (in `patch_pygeoapi.py`): +- Empty results returned `"bbox": [NaN, NaN, NaN, NaN]`, which isn't valid JSON. The + Parquet provider now leaves out both the collection bbox and the per-feature bboxes. +- `numberMatched` was an estimate except on the last page. Items pages now carry the true + total, so clients can fetch pages in parallel. +- Not fixed here (separate ticket): WQP `source` labels differ between products. See T6. + +Until T5 is done, **each pipeline run that changes data needs a manual pygeoapi build** +(`orchestration/pygeoapi/cloudbuild.yaml`), started after the product runs finish. +Whatever account runs it needs `storage.objectViewer` on the products bucket for the +bake, as well as what builds need today. + +New finding — **memory at 2 GiB with 4 workers.** Under sustained 10,000-row page load the +local container's cgroup peak reached its 2 GiB cap (anon 1.88 GiB; workers 470–650 MB +each). This is true of the live config too. Cloud Run has no swap, so this would be an OOM +restart. Same load: `ARROW_DEFAULT_MEMORY_POOL=system` peaked at 1,982 MiB (no real help); +`WSGI_WORKERS=2` peaked at 1,781 MiB with the same timings. **Decided: `WSGI_WORKERS=2`**, +set in `cloudbuild.yaml` `--set-env-vars`. Watch Cloud Run memory after deploy; raise to +`--memory=4Gi` if it still approaches the cap. + +## Goal + +The water-level hydrograph client must be at least as fast against pygeoapi as it was +against GeoServer, and paging must be correct. Every layer should get faster along the way. + +## Where things stand + +- `die-pygeoapi` is live on Cloud Run: 31 collections, 2 GiB, dedicated service account, + reading products from `gs://dataservices-die-products/products/{id}/`. +- The live `nm_waterlevels_timeseries` is served from a **manually uploaded** + `latest.parquet`. It has no `feature_id` column, is sorted by `datetime`, and carries + a stray `__index_level_0__` column. The other 30 layers are GeoJSON through OGR. +- The DIE-363 branch makes the pipeline write `latest.parquet` for every product and makes + the image build choose Parquet for each product when the file is usable (point geometry, + bbox covering column, id column, and a tz-aware `datetime` for timeseries). Locally, + 28 of 31 layers move to Parquet. The 3 polygon layers stay on GeoJSON because the Parquet + provider's `bbox` is containment-only. + +## Measurements to beat + +The hydrograph client requests one well (filter by well id) in 5,000-row pages. Timed on +the longest well (id 9033, 48,809 readings): + +| | GeoServer | pygeoapi (live, today) | +|---|---|---| +| Time per 5,000-row page | ~1.0 s | ~2.5–3 s | +| Rows returned per page | 5,000 | ~9,600 (pages overlap) | +| Uncompressed size per page | 2.8 MB | 6.9 MB | +| Gzip | yes | no | +| `numberMatched` | true total (48,809) | lookahead (offset + returned + 1) | +| Whole series | ~10 s (10 pages) | ~26 s (10 pages) | + +GeoServer is **not** backed by PostGIS. The `/geoserver` asset uploads a +GeoPackage that GeoServer stores on its own disk. Its advantages are local disk, a +long-lived datastore, and SQL paging. + +Local results (DIE-363 branch, files on local disk, 2 GiB): every layer answers items, +`bbox`, and item-by-id in 12–240 ms. Peak memory is 1,733 MiB. On the live service the same +Parquet reads cost ~1.2 s per request, because pygeoapi builds a new provider on every +request (`load_plugin` in `pygeoapi/api/itemtypes.py`), and opening a `gs://` dataset +takes several GCS round trips. + +## Root causes (verified in the pygeoapi 0.25.dev0 source in the image) + +1. **Per-request GCS overhead (~1.2 s floor).** The provider is constructed per request. + Opening `pyarrow.dataset.dataset("gs://...")`, reading the footer, and `get_fields()` + all go over the network, even for `limit=10` with no filter. +2. **Paging bug in the Parquet provider.** In + `pygeoapi/provider/parquet.py::_response_feature_collection`: + ```python + for batch in batches: + read += batch.num_rows + if read > limit: + batches_list.append(batch.slice(0, limit + 1)) # ignores rows already collected + break + ``` + When a filter spreads matches across several batches, a page holds `limit + 1` rows + *plus* every earlier batch. That's the ~9,600 rows for `limit=5000`, a constant + +~4,600 at `limit=10000`, and the overlap between pages. Fix: slice + `limit + 1 - (read - batch.num_rows)`. +3. **`numberMatched` is lookahead by design.** It's `offset + len(rows)` with one extra row + fetched to decide on a `next` link. Fixing (2) does **not** make it a true total. + (Since patched to count the filtered rows; see Progress.) +4. **No gzip.** Cloud Run doesn't compress, and `server.gzip` isn't set. +5. **Well lookups scan the whole file.** The timeseries Parquet is sorted by `datetime`, + so an `id=` filter can't skip row groups. +6. **Bigger rows.** Every Parquet column becomes a property, including the `bbox` covering + column. geopandas `__geo_interface__` also adds a bbox to every feature. The manual + file also has `__index_level_0__`, which the pipeline converter doesn't write + (`index=False`). + +## Decisions + +- **Bake the product files into the image at build time, and rebuild after each pipeline + run.** The pipeline is moving to a weekly cadence, so data only changes at rebuilds. + This removes cause (1) entirely. + - Rejected: a Cloud Storage FUSE volume mount (the speed gain depends on caching and + is unmeasured), copying at container startup (unnecessary with weekly data), and + PostGIS (new always-on infrastructure; reverses SPEC's "no database"). +- **Keep the build-time Parquet check** (per product, falling back to GeoJSON). Now it + also decides which file gets baked. +- **Patch the Parquet provider in the image** rather than fork pygeoapi, and report the + bug upstream. +- **Out of scope here:** EDR (`/locations/{id}` CoverageJSON; see "Later"), PostGIS. + +## Tasks + +Do them in order. Each is independently shippable. Re-run the measurement protocol +(below) after T3 and after T5. + +### T1 — Gzip responses + +- `orchestration/pygeoapi/config.yml.j2`: add `gzip: true` under `server:`. +- Accept: `curl -sI -H 'Accept-Encoding: gzip' .../items?limit=5000` shows + `Content-Encoding: gzip`, and the compressed size is a fraction of the raw size. + Clients that don't send `Accept-Encoding` still get plain JSON. + +### T2 — Pin the base image, patch the Parquet provider + +- `orchestration/pygeoapi/Dockerfile`: replace `geopython/pygeoapi:latest` with a pinned + digest (`geopython/pygeoapi@sha256:...`). Record the pygeoapi version it contains in a + comment. Current `:latest` is 0.25.dev0. +- Add a build step that patches the paging line (root cause 2). The patch **must fail the + build** if the expected source line isn't found, so a base-image bump can't drop it + silently. A small Python script is safer than `sed`: assert the old text occurs + exactly once, then replace it. +- Optional, same patch step: coerce naive `datetime_` bounds to UTC before comparing, so + date-only queries (`2020-01-01/2020-12-31`) stop returning 500 against the tz-aware + column. +- Open an upstream pygeoapi issue or PR for the paging bug, and link it in the Dockerfile + comment. +- Accept: + - `items?id=9033&limit=5000&offset=N` returns exactly 5,000 rows for every page but the + last, and the pages don't overlap. + - The union of all pages has 48,809 distinct `feature_id`s. + - The date-only `datetime` query returns 200 (if the optional fix is done). + +### T3 — Sort timeseries Parquet by well, then time + +- `backend/persisters/geodataframe.py::collection_to_geoparquet` (currently + `sort_values(PARQUET_TIME_FIELD)`): for timeseries collections, sort by + `("id", "source", "datetime")` (see Progress for why `id` first). Other products need + no particular order. + - Trade-off: `datetime` queries across all wells lose row-group skipping. They stay + correct, just slower. The hydrograph is the priority. + - Consider a smaller `row_group_size` (e.g. 50,000) so one well spans few row groups. + Measure first. +- Update `tests/test_persisters/test_geodataframe.py::test_datetime_is_utc_and_sorted`. +- Takes effect only when the pipeline rewrites the timeseries Parquet. Existing files + won't be rewritten on their own: `_sync_parquet` skips when the GeoJSON content hash + matches, and a converter change doesn't change that hash. Do T6's + `PARQUET_FORMAT_VERSION` first, or delete the object once. + +### T4 — Bake product files into the image + +- Build flow (replaces the GCS paths in the config): + 1. With `--network=cloudbuild` (already set), a new build step downloads, for each + product in `products.yaml`, `latest.parquet` when it passes + `generate_config.parquet_id_field`, otherwise `latest.geojson`, into + `/data/products/{id}/`. `generate_config.py` can do the download, e.g. with + `--bake-to /data/products`. + - Permission or network errors fail the build (as the check does today). + - A product with neither usable file fails the build. Better to keep the old revision + serving than to ship an image that won't start. + 2. Render the config against the baked copies (`--source-root /data/products`). +- Local builds without the build arg keep today's GCS-path config, so a local + `docker build` still works. +- Image size: ~81 MB of Parquet (all 28 point layers) plus ~15 MB of GeoJSON for the 3 + polygon layers. +- Runtime no longer reads GCS for products. `CPL_MACHINE_IS_GCE` and the service + account's bucket read become unnecessary but harmless. Remove them in a later cleanup. +- Update `SPEC.md` §6 and the invariants: sources are baked copies taken from GCS at build + time. Drop "OGR MUST use /vsigs/" and "Parquet MUST use gs://". +- Accept: + - The live floor for `limit=10` drops from ~1.2 s to local-disk levels (target < 100 ms + server time). + - The build log lists each product's baked file and the reason. + - All 31 collections return 200. + +### T5 — Rebuild after a pipeline run (deferred; manual for now) + +Not being done yet. Until it is, rebuild by hand after product runs publish new data +(SPEC §6.6). Design notes for when it's picked up: + +- Create a Cloud Build **trigger** on the GitHub repo, `main` branch, manual invocation, + using `orchestration/pygeoapi/cloudbuild.yaml`. Triggers fill `COMMIT_SHA` + automatically. +- Start it when a pipeline run has published new data. Recommended: a Dagster + job/sensor that runs after the scheduled product runs complete. If any product's + materialization has `parquet_status: written` or `skipped_unchanged: false`, it calls + the Cloud Build `projects.triggers.run` API. + - One rebuild per run, not per product. Debounce, and never build mid-run, so the image + never bakes a half-updated set. + - The Dagster+ credentials (`GCP_SERVICE_ACCOUNT_KEY`) need permission to run the + trigger (`roles/cloudbuild.builds.editor`, or a narrower custom role). + - Simpler fallback: Cloud Scheduler runs the trigger on a fixed cadence after the + pipeline window. +- Schedules today are **monthly and staggered** (`products.yaml`: days 1, 15, 22, 23, 25 + at 06:00 UTC). Moving to weekly means updating those `schedule:` values. With the + current staggering, the trigger fires up to 5 times a month. +- Accept: a pipeline run that changes a product results in exactly one build, and the new + revision serves the new data with no manual step. + +### T6 — Pipeline correctness fixes + +- **Unique GeoJSON feature ids.** `backend/persisters/ogc_features.py:565` builds + timeseries ids as `f"{source}:{site_id}:{date}"`, which repeats for multiple readings on + one day (499,711 distinct ids in 882,280 features). Use the full `datetime` (or append a + sequence number). The Parquet `feature_id` already de-duplicates with `:2`, `:3` + suffixes. Once the GeoJSON ids are unique, those suffixes stop appearing. GeoServer and + GeoJSON consumers see the id change. +- **Well id isn't unique across sources** (e.g. `10493` is ST2/BernCo). The hydrograph + client filters by `id` alone. Either filter by `source` too, or add a `site_key` + (`source:id`) column and document it. +- **Force a Parquet rewrite when the converter changes.** `_sync_parquet` + (`orchestration/resources/gcs.py`) skips when the content hash matches, so converter + changes (T3) don't reach existing files. Include a `PARQUET_FORMAT_VERSION` in the hash + stored on `latest.parquet`. +- **Data quality:** one timeseries observation is dated 2316-03-02. Find the source + record, and decide whether to reject future dates at transform time. +- **Dagster+ memory:** converting the 733 MB timeseries peaked at ~3.6 GB RSS locally + (JSON load + GeoDataFrame). Watch the first run. If it's too high, build the Parquet + from the in-memory records instead of re-parsing the GeoJSON. +- **WQP source labels differ between products** (separate ticket). In + `backend/connectors/wqp/transformer.py`, WQP site records are labeled `WQP/{provider}` + (`WQP/NWIS`, `WQP/STORET`), but observation records (water levels, analytes) are + labeled plain `WQP`. Products built from observations say `WQP`: the summaries, the + timeseries, major chemistry, and MCL exceedance. Products built from sites say + `WQP/NWIS`. The client's fallback copes, but for those wells the summary shows the + USGS-NWIS feed's row. The recommended fix is to label observations `WQP/{provider}` too. + That changes public `source` values in every WQP-backed product, on GeoServer as well. +- **CI:** `uv sync --extra dev` has no pyarrow, jinja2, or dagster, so the new Parquet, + config, and GCS-sync tests skip in CI. Add `--extra parquet` and the needed packages to + CI, or run those tests in `orchestration-ci.yml`. + +### T7 — Trim response rows (optional, after T1–T5) + +- Extend the T2 patch to drop the `bbox` covering column from properties. (Per-feature + `bbox` is already gone.) Gzip already recovers most of these bytes, + so do this only if measurements still show a gap. +- Don't drop the Z coordinate for size reasons; that's a data decision. + +## Measurement protocol + +Run it against the live service after T3 (with a build of that code) and after T5. + +1. **Timing:** time each page of `items?id=9033&limit=5000&offset=0,5000,...` until a short + page, with `Accept-Encoding: gzip`. Record the time per page, the rows per page, the + bytes on the wire (`%{size_download}`), and the total. +2. **Correctness:** the union of the pages is exactly 48,809 distinct `feature_id`s. +3. **Baseline:** `limit=10` with no filter on 3 layers (the floor). +4. **All layers:** items, `bbox`, and item-by-id for all 31 collections return 200. The + polygon layers return 5 for the standard bbox (`-107,32,-105,34`). +5. **Targets:** whole series ≤ 10 s (GeoServer parity); floor < 100 ms server time; no + OOM at 2 GiB. + +## Later + +- **EDR for the timeseries.** A `BaseEDRProvider` subclass serving + `/collections/nm_waterlevels_timeseries/locations/{source:id}` as a CoverageJSON + `PointSeries`. It sits beside the Features provider on the same collection, and returns + one request per well with no paging. Model it on OcotilloAPI's `core/edr_provider.py` + (PostgreSQL-backed). The DIE version would read the baked Parquet with a pyarrow filter + (fast after T3). Notes from that code: implement `instances()`/`instance()` exactly as + pygeoapi calls them, and return fresh dicts from `get_fields()`. Per-observation fields + such as `approval_status` and `qualifier` need to be separate parameters, or stay on + the Features endpoint. Needs NewWeaver to read CoverageJSON. +- **PostGIS**, if a shared Postgres becomes available. It makes the Ocotillo EDR provider + reusable nearly as-is, and gives SQL paging with true counts. +- **Upstream pygeoapi:** the OpenAPI contact address bug (`gen_contact` overwrites + `administrativeArea` with postal code and country). diff --git a/orchestration/assets/products.py b/orchestration/assets/products.py index cf2d317..ffd2787 100644 --- a/orchestration/assets/products.py +++ b/orchestration/assets/products.py @@ -547,6 +547,7 @@ def _combine_asset( # dated_uri is None when the upload was skipped as unchanged. if info.get("dated_uri"): metadata["dated_uri"] = dg.MetadataValue.url(info["dated_uri"]) + metadata.update(_parquet_metadata(info)) # Change-recency signal for tuning run frequency: how long the data has # been static. A large, growing value means the schedule can be relaxed # (e.g. daily -> monthly); 0 means it changed this run. @@ -637,6 +638,7 @@ def _combine_asset( ) if info.get("dated_uri"): metadata["dated_uri"] = dg.MetadataValue.url(info["dated_uri"]) + metadata.update(_parquet_metadata(info)) yield dg.Output(None, metadata=metadata) yield dg.AssetCheckResult( @@ -730,6 +732,7 @@ def _combine_asset( ) if info.get("dated_uri"): metadata["dated_uri"] = dg.MetadataValue.url(info["dated_uri"]) + metadata.update(_parquet_metadata(info)) yield dg.Output(None, metadata=metadata) yield dg.AssetCheckResult( @@ -746,6 +749,19 @@ def _combine_asset( return _combine_asset +def _parquet_metadata(info: dict) -> dict: + """Metadata for latest.parquet. A failed conversion doesn't fail the + product, so it's surfaced here.""" + meta: dict = {} + if info.get("parquet_status"): + meta["parquet_status"] = dg.MetadataValue.text(info["parquet_status"]) + if info.get("parquet_status") != "failed" and info.get("parquet_uri"): + meta["parquet_uri"] = dg.MetadataValue.url(info["parquet_uri"]) + if info.get("parquet_error"): + meta["parquet_error"] = dg.MetadataValue.text(info["parquet_error"]) + return meta + + def _num_opt(value): """None passes through; otherwise coerce to float (products.yaml overrides for the correlation thresholds are optional).""" diff --git a/orchestration/pygeoapi/Dockerfile b/orchestration/pygeoapi/Dockerfile index 5c6a088..4127b4a 100644 --- a/orchestration/pygeoapi/Dockerfile +++ b/orchestration/pygeoapi/Dockerfile @@ -1,7 +1,8 @@ -FROM geopython/pygeoapi:latest +# Pinned for patch_pygeoapi.py: pygeoapi 0.25.dev0 (:latest on 2026-09-24). +FROM geopython/pygeoapi@sha256:f3dd50a56f870d80df67416b6b07bf4d92aefa78b7e6f697564288f1cdf1df3f -# pygeoapi base image includes GDAL with /vsigs/ GCS support. -# Auth uses Application Default Credentials — no key file needed on Cloud Run. +# Unbaked images read GCS directly (GDAL /vsigs/, Application Default +# Credentials). WORKDIR /pygeoapi @@ -11,6 +12,10 @@ WORKDIR /pygeoapi RUN /venv/bin/python3 -m pip install --no-cache-dir --quiet \ pyarrow==25.0.1 geopandas==1.1.4 +# Fix Parquet provider bugs; fails the build if the patched source changed. +COPY pygeoapi/patch_pygeoapi.py /tmp/patch_pygeoapi.py +RUN /venv/bin/python3 /tmp/patch_pygeoapi.py + # Copy generation inputs # Build context is orchestration/ (see cloudbuild.yaml) COPY pygeoapi/config.yml.j2 /tmp/config.yml.j2 @@ -20,10 +25,19 @@ COPY config/products.yaml /tmp/products.yaml # Bake config into image at build time. # §V: config generated from products.yaml, not hand-edited. # jinja2 and PyYAML are pygeoapi dependencies, already in /venv. +# BAKE_PRODUCTS=1 (cloudbuild.yaml): copy each product's file from GCS into +# /data/products and serve it locally; needs GCS access at build time. Use +# --no-cache locally to re-download. +# CHECK_PARQUET=1: choose Parquet vs GeoJSON the same way, but read GCS. +# Neither: GeoJSON from GCS. +ARG BAKE_PRODUCTS=0 +ARG CHECK_PARQUET=0 RUN /venv/bin/python3 /tmp/generate_config.py \ --products /tmp/products.yaml \ --template /tmp/config.yml.j2 \ - --output /pygeoapi/local.config.yml + --output /pygeoapi/local.config.yml \ + $([ "$CHECK_PARQUET" = "1" ] && echo --check-parquet) \ + $([ "$BAKE_PRODUCTS" = "1" ] && echo --bake-to /data/products) EXPOSE 80 diff --git a/orchestration/pygeoapi/cloudbuild.yaml b/orchestration/pygeoapi/cloudbuild.yaml index b7bcb32..fc65218 100644 --- a/orchestration/pygeoapi/cloudbuild.yaml +++ b/orchestration/pygeoapi/cloudbuild.yaml @@ -2,6 +2,10 @@ steps: - name: 'gcr.io/cloud-builders/docker' args: - 'build' + # Bakes products from GCS into the image; the cloudbuild network supplies + # credentials (needs storage.objectViewer on the products bucket). + - '--network=cloudbuild' + - '--build-arg=BAKE_PRODUCTS=1' - '-t' - '$_REGION-docker.pkg.dev/$PROJECT_ID/die/die-pygeoapi:$COMMIT_SHA' - '-t' @@ -28,10 +32,11 @@ steps: - '--allow-unauthenticated' - '--port=80' - '--service-account=die-pygeoapi@$PROJECT_ID.iam.gserviceaccount.com' # needs storage.objectViewer on the products bucket - - '--memory=2Gi' # Parquet provider peaked at 1.2-2.0 GB across 4 workers locally + - '--memory=2Gi' # CPL_MACHINE_IS_GCE: GDAL can't detect Cloud Run as GCE, so /vsigs/ skips the # metadata-server credentials and every GeoJSON collection fails without it. - - '--set-env-vars=PYGEOAPI_SERVER_URL=https://die-pygeoapi-$PROJECT_NUMBER.$_REGION.run.app,CPL_MACHINE_IS_GCE=YES' + # WSGI_WORKERS=2: the default 4 hit the 2 GiB cap under load locally. + - '--set-env-vars=PYGEOAPI_SERVER_URL=https://die-pygeoapi-$PROJECT_NUMBER.$_REGION.run.app,CPL_MACHINE_IS_GCE=YES,WSGI_WORKERS=2' substitutions: _REGION: us-central1 diff --git a/orchestration/pygeoapi/config.yml.j2 b/orchestration/pygeoapi/config.yml.j2 index ea00c1a..a1b85c6 100644 --- a/orchestration/pygeoapi/config.yml.j2 +++ b/orchestration/pygeoapi/config.yml.j2 @@ -8,6 +8,7 @@ server: language: en-US cors: true pretty_print: false + gzip: true # Cloud Run doesn't compress limits: default_items: 500 max_items: 10000 # NewWeaver requests limit=10000 @@ -73,22 +74,22 @@ resources: end: null {% endif %} providers: -{% if product.output_type == 'ogc_timeseries' %} - # GeoParquet via the Parquet provider: the OGR provider ignores datetime - # filters and parses the whole GeoJSON on every request. The file must be - # GeoParquet with a bbox covering column and a tz-aware UTC datetime. +{% if product.id in parquet %} + # Chosen by generate_config.py when latest.parquet is usable. - type: feature name: Parquet data: - source: gs://{{ gcs_bucket }}/products/{{ product.id }}/latest.parquet - id_field: id + source: {{ parquet_root }}/{{ product.id }}/latest.parquet + id_field: {{ parquet[product.id] }} +{% if product.output_type == 'ogc_timeseries' %} time_field: datetime +{% endif %} {% else %} - type: feature name: OGR data: source_type: GeoJSON - source: /vsigs/{{ gcs_bucket }}/products/{{ product.id }}/latest.geojson + source: {{ geojson_root }}/{{ product.id }}/latest.geojson source_options: GDAL_HTTP_UNSAFESSL: "NO" gdal_ogr_options: diff --git a/orchestration/pygeoapi/generate_config.py b/orchestration/pygeoapi/generate_config.py index 4610200..f4d2215 100644 --- a/orchestration/pygeoapi/generate_config.py +++ b/orchestration/pygeoapi/generate_config.py @@ -2,24 +2,136 @@ Generate pygeoapi config.yml from products.yaml + Jinja2 template. §V: pygeoapi config MUST be generated from products.yaml — never hand-edited. -§V: pygeoapi OGR provider MUST use /vsigs/ path (GCS). Timeseries products use - the Parquet provider and read gs://.../latest.parquet instead. + +--check-parquet: serve each product from latest.parquet when usable, else +latest.geojson. --bake-to DIR (implies --check-parquet): also copy the chosen +file into DIR and serve it from there. --source-root: read a local directory +instead of GCS (testing only). Usage: python generate_config.py \ --products ../config/products.yaml \ --template config.yml.j2 \ - --output /pygeoapi/local.config.yml + --output /pygeoapi/local.config.yml \ + [--check-parquet | --bake-to /data/products] [--source-root DIR] """ import argparse +import json from pathlib import Path +from typing import Optional import yaml from jinja2 import Environment, FileSystemLoader +PARQUET_ID_FIELDS = ("feature_id", "id") # first found wins; feature_id is unique +TIME_FIELD = "datetime" +# The Parquet provider's bbox filter is containment-only, which drops polygons +# crossing the box edge, so only point layers use it. +POINT_TYPES = {"Point", "Point Z"} + + +def parquet_id_field(uri: str, needs_time: bool) -> tuple[Optional[str], str]: + """Check the Parquet footer at *uri* (gs:// or local). Returns + ``(id_field, "ok")`` if usable, else ``(None, reason)``. + + Usable: GeoParquet metadata with a bbox covering column, point geometry, + an id column, and (if *needs_time*) a tz-aware ``datetime``. A missing or + unreadable file returns None; other errors (e.g. permissions) raise, so + the build fails instead of silently serving GeoJSON.""" + import pyarrow as pa + import pyarrow.fs as pafs + import pyarrow.parquet as pq + + fs, path = pafs.FileSystem.from_uri(uri) + if fs.get_file_info(path).type == pafs.FileType.NotFound: + return None, "no latest.parquet" + try: + schema = pq.read_schema(path, filesystem=fs) + except pa.ArrowInvalid as exc: + return None, f"unreadable Parquet ({exc})" + + geo = json.loads((schema.metadata or {}).get(b"geo", b"null")) + if not geo: + return None, "no GeoParquet metadata" + geom = geo.get("columns", {}).get(geo.get("primary_column"), {}) + if not geom.get("covering", {}).get("bbox"): + return None, "no bbox covering column" + types = set(geom.get("geometry_types") or []) + if not types or not types <= POINT_TYPES: + return None, f"non-point geometry {sorted(types)} (Parquet bbox is containment-only)" + if needs_time: + if TIME_FIELD not in schema.names: + return None, f"no {TIME_FIELD} column" + t = schema.field(TIME_FIELD).type + if not (pa.types.is_timestamp(t) and t.tz): + return None, f"{TIME_FIELD} is {t}, not a tz-aware timestamp" + for field in PARQUET_ID_FIELDS: + if field in schema.names: + return field, "ok" + return None, "no id column" + + +def bake_file(uri: str, dest: Path) -> int: + """Copy *uri* (gs:// or a local path) to *dest*; returns its size in bytes.""" + import pyarrow.fs as pafs + + fs, path = pafs.FileSystem.from_uri(uri) + dest.parent.mkdir(parents=True, exist_ok=True) + pafs.copy_files( + path, str(dest), source_filesystem=fs, destination_filesystem=pafs.LocalFileSystem() + ) + return dest.stat().st_size + + +def _exists(uri: str) -> bool: + import pyarrow.fs as pafs + + fs, path = pafs.FileSystem.from_uri(uri) + return fs.get_file_info(path).type != pafs.FileType.NotFound -def generate(products_path: Path, template_path: Path, output_path: Path) -> None: + +def generate( + products_path: Path, + template_path: Path, + output_path: Path, + check_parquet: bool = False, + source_root: Optional[str] = None, + bake_to: Optional[Path] = None, +) -> None: products_config = yaml.safe_load(products_path.read_text()) + bucket = products_config["gcs_bucket"] + store_root = source_root or f"gs://{bucket}/products" + if bake_to: + check_parquet = True + geojson_root = parquet_root = str(bake_to) + else: + geojson_root = source_root or f"/vsigs/{bucket}/products" + parquet_root = store_root + + # product id -> Parquet id_field, for products served from GeoParquet + parquet: dict[str, str] = {} + for product in products_config["products"]: + pid = product["id"] + if check_parquet: + id_field, reason = parquet_id_field( + f"{store_root}/{pid}/latest.parquet", + needs_time=product.get("output_type") == "ogc_timeseries", + ) + else: + id_field, reason = None, "Parquet check disabled" + if id_field: + parquet[pid] = id_field + name, line = "latest.parquet", f"Parquet (id_field={id_field})" + else: + name, line = "latest.geojson", f"GeoJSON ({reason})" + if bake_to: + # Fail the build so the previous revision keeps serving. + if not _exists(f"{store_root}/{pid}/{name}"): + raise SystemExit(f"{pid}: no {name} to bake ({line})") + size = bake_file(f"{store_root}/{pid}/{name}", Path(bake_to) / pid / name) + line += f", baked {size / 1e6:.1f} MB" + print(f" {pid}: {line}") + env = Environment( loader=FileSystemLoader(str(template_path.parent)), keep_trailing_newline=True, @@ -27,24 +139,28 @@ def generate(products_path: Path, template_path: Path, output_path: Path) -> Non tmpl = env.get_template(template_path.name) rendered = tmpl.render( products=products_config["products"], - gcs_bucket=products_config["gcs_bucket"], + parquet=parquet, + geojson_root=geojson_root, + parquet_root=parquet_root, ) - # Sanity check: every product must produce a GCS source entry + # Sanity check: every product must produce a source entry for its provider for product in products_config["products"]: pid = product["id"] - bucket = products_config["gcs_bucket"] - if product.get("output_type") == "ogc_timeseries": - expected = f"gs://{bucket}/products/{pid}/latest.parquet" + if pid in parquet: + expected = f"{parquet_root}/{pid}/latest.parquet" else: - expected = f"/vsigs/{bucket}/products/{pid}/latest.geojson" + expected = f"{geojson_root}/{pid}/latest.geojson" assert expected in rendered, ( f"§V violated: provider for '{pid}' missing {expected} in generated config" ) output_path.parent.mkdir(parents=True, exist_ok=True) output_path.write_text(rendered) - print(f"Generated {output_path} ({len(products_config['products'])} collections)") + print( + f"Generated {output_path} ({len(products_config['products'])} collections, " + f"{len(parquet)} from GeoParquet)" + ) if __name__ == "__main__": @@ -52,5 +168,26 @@ def generate(products_path: Path, template_path: Path, output_path: Path) -> Non parser.add_argument("--products", required=True, type=Path) parser.add_argument("--template", required=True, type=Path) parser.add_argument("--output", required=True, type=Path) + parser.add_argument( + "--check-parquet", + action="store_true", + help="serve each product from latest.parquet when a usable one exists", + ) + parser.add_argument( + "--bake-to", + type=Path, + help="copy each product's chosen file here and serve it from local disk", + ) + parser.add_argument( + "--source-root", + help="local testing only: read products from this directory instead of GCS", + ) args = parser.parse_args() - generate(args.products, args.template, args.output) + generate( + args.products, + args.template, + args.output, + args.check_parquet, + args.source_root, + args.bake_to, + ) diff --git a/orchestration/pygeoapi/patch_pygeoapi.py b/orchestration/pygeoapi/patch_pygeoapi.py new file mode 100644 index 0000000..be4d888 --- /dev/null +++ b/orchestration/pygeoapi/patch_pygeoapi.py @@ -0,0 +1,83 @@ +""" +Patch pygeoapi's Parquet provider at build time. Each target must occur exactly +once, or the build fails. Delete a fix once it lands upstream. + +1. Paging: the batch crossing `limit` was sliced to limit + 1 rows, ignoring + rows already read, so pages were too long and overlapped. +2. Date-only `datetime` bounds were naive and raised against the tz-aware + column (500). Treat them as UTC; time_field is only set for tz-aware files. +3. bbox: `__geo_interface__` added bboxes, [NaN, NaN, NaN, NaN] (invalid JSON) + when empty. Omit them. +4. numberMatched: report the filtered row count, not offset + rows returned. + +Usage: + python patch_pygeoapi.py [path/to/pygeoapi/provider/parquet.py] +""" +import sys +from pathlib import Path + +FIXES = [ + ( + "paging: count rows already read", + "batches_list.append(batch.slice(0, limit + 1))", + "batches_list.append(\n" + " batch.slice(0, limit + 1 - (read - batch.num_rows)))", + ), + ("datetime: naive begin is UTC", "begin = isoparse(begin)", "begin = _utc(isoparse(begin))"), + ("datetime: naive end is UTC", "end = isoparse(end)", "end = _utc(isoparse(end))"), + ( + "datetime: naive instant is UTC", + "target_time = isoparse(datetime_)", + "target_time = _utc(isoparse(datetime_))", + ), + ( + "bbox: items response", + "result = gdf.__geo_interface__", + "result = gdf.to_geo_dict(na='null', show_bbox=False, drop_id=False)", + ), + ( + "bbox: single item", + "return gdf.__geo_interface__['features'][0]", + "return gdf.to_geo_dict(\n" + " na='null', show_bbox=False, drop_id=False)['features'][0]", + ), + ( + "numberMatched: true total", + "number_matched = offset + len(rp)", + "number_matched = scanner.count_rows()", + ), +] + +UTC_HELPER = ''' + +def _utc(value): + """Added by DIE's patch_pygeoapi.py: naive datetimes are UTC.""" + from datetime import timezone + return value if value.tzinfo else value.replace(tzinfo=timezone.utc) +''' + + +def patch(source: str) -> str: + for name, old, new in FIXES: + count = source.count(old) + if count != 1: + raise SystemExit( + f"patch '{name}': expected 1 occurrence of {old!r}, found {count}. " + "pygeoapi changed; re-check the fix against the new source." + ) + source = source.replace(old, new) + return source + UTC_HELPER + + +def default_path() -> Path: + import importlib.util + + # find_spec locates the package without importing the provider module. + spec = importlib.util.find_spec("pygeoapi") + return Path(spec.origin).parent / "provider" / "parquet.py" + + +if __name__ == "__main__": + path = Path(sys.argv[1]) if len(sys.argv) > 1 else default_path() + path.write_text(patch(path.read_text())) + print(f"Patched {path} ({len(FIXES)} fixes)") diff --git a/orchestration/resources/gcs.py b/orchestration/resources/gcs.py index a118ac0..bdd145a 100644 --- a/orchestration/resources/gcs.py +++ b/orchestration/resources/gcs.py @@ -16,6 +16,7 @@ _CONTENT_HASH_KEY = "content_hash" _LAST_CHANGED_KEY = "last_changed" # YYYY-MM-DD the content last actually changed +_PARQUET_FORMAT_KEY = "parquet_format" # converter version that wrote latest.parquet def _days_between(start: str, end: str) -> Optional[int]: @@ -28,13 +29,11 @@ def _days_between(start: str, end: str) -> Optional[int]: return None -def _content_hash(local_path: str) -> str: - """Stable SHA-256 of a product's GeoJSON content, ignoring the volatile +def _content_hash(data: dict) -> str: + """Stable SHA-256 of a product's parsed GeoJSON *data*, ignoring the volatile `timeStamp` so that re-running with unchanged data yields the same hash.""" - with open(local_path, encoding="utf-8") as f: - data = json.load(f) - data.pop("timeStamp", None) - payload = json.dumps(data, sort_keys=True, default=str).encode("utf-8") + stable = {k: v for k, v in data.items() if k != "timeStamp"} + payload = json.dumps(stable, sort_keys=True, default=str).encode("utf-8") return hashlib.sha256(payload).hexdigest() @@ -53,7 +52,8 @@ def _storage_client(): class GCSResource(dg.ConfigurableResource): """ - Upload OGC Feature Collection GeoJSON files to GCS. + Upload OGC Feature Collection GeoJSON files to GCS, each with a GeoParquet + copy (latest.parquet) for pygeoapi's Parquet provider. §V: latest.geojson MUST be overwritten atomically (copy from dated object, never direct overwrite of in-flight file). @@ -95,6 +95,8 @@ def upload_product( Dedup compares a content hash (ignoring the volatile timeStamp) against the hash stored on the current latest.geojson's metadata. + Either way, latest.parquet is then synced (:meth:`_sync_parquet`). + Returns dict with: dated_uri: gs://bucket/products/{product_id}/{date}.geojson (None if skipped) latest_uri: gs://bucket/products/{product_id}/latest.geojson @@ -102,6 +104,9 @@ def upload_product( file_size_bytes: int run_date: str skipped: bool — True when content matched and nothing was written + parquet_uri: gs://bucket/products/{product_id}/latest.parquet + parquet_status: "written" | "unchanged" | "failed" + parquet_error: str — only when parquet_status is "failed" """ if run_date is None: run_date = datetime.now(timezone.utc).strftime("%Y-%m-%d") @@ -113,10 +118,9 @@ def upload_product( latest_key = f"{self.products_prefix}/{product_id}/latest.geojson" file_size = Path(local_path).stat().st_size - new_hash = _content_hash(local_path) - with open(local_path, encoding="utf-8") as f: data = json.load(f) + new_hash = _content_hash(data) feature_count = data.get("numberReturned", len(data.get("features", []))) latest_uri = f"gs://{self.bucket_name}/{latest_key}" @@ -140,6 +144,7 @@ def upload_product( "skipped": True, "last_changed": last_changed, "days_since_last_change": _days_between(last_changed, run_date), + **self._sync_parquet(bucket, product_id, data, new_hash, local_path), } # Content changed (or first upload) — last_changed is now. @@ -161,8 +166,49 @@ def upload_product( "skipped": False, "last_changed": run_date, "days_since_last_change": 0, + **self._sync_parquet(bucket, product_id, data, new_hash, local_path), } + def _sync_parquet( + self, bucket, product_id: str, data: dict, content_hash: str, local_path: str + ) -> dict: + """Write latest.parquet from *data* unless it already carries + *content_hash* and the current PARQUET_FORMAT_VERSION. A failed + conversion doesn't fail the product (GeoJSON is the fallback), but + deletes any stale latest.parquet.""" + key = f"{self.products_prefix}/{product_id}/latest.parquet" + info: dict = {"parquet_uri": f"gs://{self.bucket_name}/{key}"} + from backend.persisters.geodataframe import ( + PARQUET_FORMAT_VERSION, + collection_to_geoparquet, + ) + + wanted = {_CONTENT_HASH_KEY: content_hash, _PARQUET_FORMAT_KEY: PARQUET_FORMAT_VERSION} + blob = bucket.blob(key) + if blob.exists(): + blob.reload() + meta = blob.metadata or {} + if all(meta.get(k) == v for k, v in wanted.items()): + return {**info, "parquet_status": "unchanged"} + + parquet_path = Path(local_path).with_name(f"{product_id}.parquet") + try: + collection_to_geoparquet(data, parquet_path) + blob.metadata = wanted + blob.upload_from_filename( + str(parquet_path), content_type="application/vnd.apache.parquet" + ) + except Exception as exc: # noqa: BLE001 — soft-fail, GeoJSON is the fallback + dg.get_dagster_logger().warning( + f"{product_id}: GeoParquet not written, pygeoapi will serve GeoJSON: {exc}" + ) + if blob.exists(): + blob.delete() + return {**info, "parquet_status": "failed", "parquet_error": str(exc)} + finally: + parquet_path.unlink(missing_ok=True) + return {**info, "parquet_status": "written"} + from dagster_gcp.gcs import GCSResource as _DagsterGCSResource # noqa: E402 diff --git a/tests/test_gcs_parquet_sync.py b/tests/test_gcs_parquet_sync.py new file mode 100644 index 0000000..c59a092 --- /dev/null +++ b/tests/test_gcs_parquet_sync.py @@ -0,0 +1,126 @@ +"""Tests for GCSResource's latest.parquet sync.""" + +import json + +import pytest + +pytest.importorskip("dagster") +pytest.importorskip("pyarrow") + +from orchestration.resources import gcs as gcs_mod # noqa: E402 + + +class FakeBlob: + def __init__(self, store, key): + self.store, self.key = store, key + self.metadata = None + + def exists(self): + return self.key in self.store + + def reload(self): + self.metadata = self.store[self.key]["metadata"] + + def upload_from_filename(self, path, content_type=None): + with open(path, "rb") as f: + self.store[self.key] = {"data": f.read(), "metadata": self.metadata} + + def delete(self): + del self.store[self.key] + + +class FakeBucket: + def __init__(self): + self.store = {} + + def blob(self, key): + return FakeBlob(self.store, key) + + def copy_blob(self, blob, bucket, new_key): + self.store[new_key] = dict(self.store[blob.key]) + + +def _resource(monkeypatch, bucket): + res = gcs_mod.GCSResource(bucket_name="bkt") + client = type("C", (), {"bucket": lambda self, name: bucket})() + monkeypatch.setattr(gcs_mod, "_storage_client", lambda: client) + return res + + +def _write_collection(tmp_path, n=2): + path = tmp_path / "collection.geojson" + features = [ + { + "type": "Feature", + "id": f"src:W{i}", + "geometry": {"type": "Point", "coordinates": [-106.5, 35.0]}, + "properties": {"id": i}, + } + for i in range(n) + ] + path.write_text(json.dumps({"type": "FeatureCollection", "features": features})) + return path + + +PARQUET = "products/p/latest.parquet" + + +def test_first_upload_writes_parquet_with_content_hash(tmp_path, monkeypatch): + bucket = FakeBucket() + info = _resource(monkeypatch, bucket).upload_product(str(_write_collection(tmp_path)), "p") + assert info["parquet_status"] == "written" + assert info["parquet_uri"] == f"gs://bkt/{PARQUET}" + geojson_hash = bucket.store["products/p/latest.geojson"]["metadata"]["content_hash"] + assert bucket.store[PARQUET]["metadata"]["content_hash"] == geojson_hash + assert bucket.store[PARQUET]["data"][:4] == b"PAR1" + + +def test_unchanged_geojson_backfills_missing_parquet_then_skips(tmp_path, monkeypatch): + bucket = FakeBucket() + res = _resource(monkeypatch, bucket) + path = str(_write_collection(tmp_path)) + res.upload_product(path, "p") + del bucket.store[PARQUET] # e.g. products published before Parquet existed + + again = res.upload_product(path, "p") + assert again["skipped"] is True + assert again["parquet_status"] == "written" + + third = res.upload_product(path, "p") + assert third["parquet_status"] == "unchanged" + + +def test_format_version_bump_rewrites_unchanged_geojson(tmp_path, monkeypatch): + import backend.persisters.geodataframe as gdf_mod + + bucket = FakeBucket() + res = _resource(monkeypatch, bucket) + path = str(_write_collection(tmp_path)) + res.upload_product(path, "p") + assert bucket.store[PARQUET]["metadata"]["parquet_format"] == gdf_mod.PARQUET_FORMAT_VERSION + assert res.upload_product(path, "p")["parquet_status"] == "unchanged" + + monkeypatch.setattr(gdf_mod, "PARQUET_FORMAT_VERSION", "next") + again = res.upload_product(path, "p") + assert again["skipped"] is True # GeoJSON unchanged + assert again["parquet_status"] == "written" + assert bucket.store[PARQUET]["metadata"]["parquet_format"] == "next" + + +def test_failed_conversion_deletes_stale_parquet(tmp_path, monkeypatch): + bucket = FakeBucket() + res = _resource(monkeypatch, bucket) + res.upload_product(str(_write_collection(tmp_path, n=2)), "p") + assert PARQUET in bucket.store + + def boom(data, out_path): + raise ValueError("bad column") + + import backend.persisters.geodataframe as gdf_mod + + monkeypatch.setattr(gdf_mod, "collection_to_geoparquet", boom) + info = res.upload_product(str(_write_collection(tmp_path, n=3)), "p") + assert info["skipped"] is False # the GeoJSON is still published + assert info["parquet_status"] == "failed" + assert "bad column" in info["parquet_error"] + assert PARQUET not in bucket.store diff --git a/tests/test_persisters/test_geodataframe.py b/tests/test_persisters/test_geodataframe.py index 9e3e8b4..be85734 100644 --- a/tests/test_persisters/test_geodataframe.py +++ b/tests/test_persisters/test_geodataframe.py @@ -19,6 +19,7 @@ dicts_to_parquet_bytes, features_to_geodataframe, geodataframe_to_features, + collection_to_geoparquet, geojson_to_geopackage, gdf_to_parquet_bytes, parquet_bytes_to_dicts, @@ -323,3 +324,98 @@ def test_converts_product_geojson(self, tmp_path): back = gpd.read_file(gpkg_path) assert len(back) == 2 assert back.geometry.iloc[0].has_z is False # flattened for GeoServer + + +class TestCollectionToGeoParquet: + """Product FeatureCollection -> latest.parquet for pygeoapi.""" + + @staticmethod + def _feature(fid, coords, **props): + geom = None if coords is None else {"type": "Point", "coordinates": coords} + return {"type": "Feature", "id": fid, "geometry": geom, "properties": props} + + def _write(self, tmp_path, features): + pytest.importorskip("pyarrow") + import pyarrow.parquet as pq + + out = tmp_path / "latest.parquet" + n = collection_to_geoparquet({"features": features}, out) + return n, pq.read_table(out) + + def test_geoparquet_metadata_and_string_columns(self, tmp_path): + n, table = self._write( + tmp_path, [self._feature("a:1", [-106.5, 35.0], name="x", value=1.5)] + ) + import pyarrow as pa + + assert n == 1 + geo = json.loads(table.schema.metadata[b"geo"]) + assert geo["version"] == "1.1.0" + assert geo["columns"]["geometry"]["covering"]["bbox"] + assert "bbox" in table.schema.names + assert not any(pa.types.is_large_string(f.type) for f in table.schema) + + def test_feature_id_is_unique_top_level_id(self, tmp_path): + _, table = self._write( + tmp_path, + [ + self._feature("src:W1:2024-01-01", [-106.5, 35.0], id=7), + self._feature("src:W1:2024-01-01", [-106.5, 35.0], id=7), + self._feature("src:W2:2024-01-01", [-106.6, 35.1], id=8), + ], + ) + rows = table.to_pydict() + assert rows["feature_id"] == [ + "src:W1:2024-01-01", + "src:W1:2024-01-01:2", + "src:W2:2024-01-01", + ] + assert rows["id"] == [7, 7, 8] # the property is kept as-is + + def test_datetime_is_utc_and_sorted(self, tmp_path): + import pyarrow as pa + + _, table = self._write( + tmp_path, + [ + self._feature("c", [-106, 35], datetime="2316-03-02T00:00:00Z"), + self._feature("a", [-106, 35], datetime="2024-01-01T06:30:00Z"), + self._feature("b", [-106, 35], datetime="2020-05-01"), + self._feature("d", [-106, 35], datetime="not a date"), + ], + ) + t = table.schema.field("datetime").type + assert pa.types.is_timestamp(t) and t.tz == "UTC" and t.unit == "ms" + values = table.to_pydict()["datetime"] + assert [v.isoformat() if v else None for v in values] == [ + "2020-05-01T00:00:00+00:00", + "2024-01-01T06:30:00+00:00", + "2316-03-02T00:00:00+00:00", # past pandas' 2262 ns limit + None, + ] + + def test_rows_sorted_by_site_then_time(self, tmp_path): + _, table = self._write( + tmp_path, + [ + self._feature("b2", [-106, 35], source="b", id="1", datetime="2021-01-01"), + self._feature("a2", [-106, 35], source="a", id="2", datetime="2020-01-01"), + self._feature("a1", [-106, 35], source="a", id="1", datetime="2022-01-01"), + self._feature("a1", [-106, 35], source="a", id="1", datetime="2019-01-01"), + ], + ) + assert table.to_pydict()["feature_id"] == ["a1:2", "a1", "b2", "a2"] + + def test_mixed_and_nested_values_become_strings(self, tmp_path): + _, table = self._write( + tmp_path, + [ + self._feature("a", [-106, 35], value=1.5, status=True, q={"k": 1}), + self._feature("b", [-106, 35], value="ND", status="Approved", q=[1, 2]), + self._feature("c", None, value=None, status=None, q=None), + ], + ) + rows = table.to_pydict() + assert rows["value"] == ["1.5", "ND", None] + assert rows["status"] == ["True", "Approved", None] + assert rows["q"] == ['{"k": 1}', "[1, 2]", None] diff --git a/tests/test_pygeoapi_config.py b/tests/test_pygeoapi_config.py new file mode 100644 index 0000000..cb4efcb --- /dev/null +++ b/tests/test_pygeoapi_config.py @@ -0,0 +1,162 @@ +"""Tests for orchestration/pygeoapi/generate_config.py.""" + +import importlib.util +from pathlib import Path + +import pytest +import yaml + +pytest.importorskip("jinja2") +pytest.importorskip("pyarrow") + +from backend.persisters.geodataframe import collection_to_geoparquet # noqa: E402 + +PYGEOAPI_DIR = Path(__file__).resolve().parents[1] / "orchestration" / "pygeoapi" +TEMPLATE = PYGEOAPI_DIR / "config.yml.j2" + + +def _load_generate_config(): + spec = importlib.util.spec_from_file_location( + "generate_config", PYGEOAPI_DIR / "generate_config.py" + ) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +gc = _load_generate_config() + + +def _point(fid, **props): + return { + "type": "Feature", + "id": fid, + "geometry": {"type": "Point", "coordinates": [-106.5, 35.0]}, + "properties": props, + } + + +def _polygon(fid): + ring = [[-107, 34], [-106, 34], [-106, 35], [-107, 35], [-107, 34]] + return { + "type": "Feature", + "id": fid, + "geometry": {"type": "Polygon", "coordinates": [ring]}, + "properties": {"name": fid}, + } + + +def _write_parquet(root: Path, pid: str, features: list) -> None: + (root / pid).mkdir(parents=True, exist_ok=True) + collection_to_geoparquet({"features": features}, root / pid / "latest.parquet") + + +def _render(tmp_path, products, check_parquet=True, source_root=None, bake_to=None) -> dict: + products_path = tmp_path / "products.yaml" + products_path.write_text(yaml.safe_dump({"gcs_bucket": "bkt", "products": products})) + out = tmp_path / "config.yml" + gc.generate(products_path, TEMPLATE, out, check_parquet, source_root, bake_to) + return yaml.safe_load(out.read_text())["resources"] + + +def _product(pid, output_type="ogc_summary"): + return {"id": pid, "title": pid, "description": pid, "output_type": output_type} + + +def test_without_check_every_product_is_geojson_on_gcs(tmp_path): + res = _render(tmp_path, [_product("a"), _product("ts", "ogc_timeseries")], False) + for pid in ("a", "ts"): + provider = res[pid]["providers"][0] + assert provider["name"] == "OGR" + assert provider["data"]["source"] == f"/vsigs/bkt/products/{pid}/latest.geojson" + + +def test_usable_point_parquet_is_served_with_feature_id(tmp_path): + root = tmp_path / "products" + _write_parquet(root, "a", [_point("x:1"), _point("x:1")]) + res = _render(tmp_path, [_product("a")], source_root=str(root)) + provider = res["a"]["providers"][0] + assert provider["name"] == "Parquet" + assert provider["data"]["source"] == f"{root}/a/latest.parquet" + assert provider["id_field"] == "feature_id" + assert "time_field" not in provider + + +def test_timeseries_parquet_gets_time_field(tmp_path): + root = tmp_path / "products" + _write_parquet(root, "ts", [_point("x:1", datetime="2024-01-01T00:00:00Z")]) + res = _render(tmp_path, [_product("ts", "ogc_timeseries")], source_root=str(root)) + provider = res["ts"]["providers"][0] + assert provider["name"] == "Parquet" + assert provider["time_field"] == "datetime" + + +@pytest.mark.parametrize( + "setup, reason", + [ + (lambda root: None, "no latest.parquet"), + (lambda root: _write_parquet(root, "a", [_polygon("county")]), "non-point"), + (lambda root: (root / "a").mkdir(parents=True) or (root / "a" / "latest.parquet").write_text("x"), + "unreadable"), + ], +) +def test_unusable_parquet_falls_back_to_geojson(tmp_path, capsys, setup, reason): + root = tmp_path / "products" + root.mkdir() + setup(root) + res = _render(tmp_path, [_product("a")], source_root=str(root)) + provider = res["a"]["providers"][0] + assert provider["name"] == "OGR" + assert provider["data"]["source"] == f"{root}/a/latest.geojson" + assert reason in capsys.readouterr().out + + +def test_timeseries_needs_tz_aware_datetime(tmp_path): + import geopandas as gpd + from shapely.geometry import Point + + root = tmp_path / "products" + (root / "ts").mkdir(parents=True) + gdf = gpd.GeoDataFrame( + {"feature_id": ["x"], "datetime": [__import__("pandas").Timestamp("2024-01-01")]}, + geometry=[Point(-106.5, 35.0)], + crs="EPSG:4326", + ) + gdf.to_parquet(root / "ts" / "latest.parquet", write_covering_bbox=True) + res = _render(tmp_path, [_product("ts", "ogc_timeseries")], source_root=str(root)) + assert res["ts"]["providers"][0]["name"] == "OGR" + + +def test_older_parquet_without_feature_id_uses_id(tmp_path): + import geopandas as gpd + from shapely.geometry import Point + + root = tmp_path / "products" + (root / "a").mkdir(parents=True) + gdf = gpd.GeoDataFrame({"id": ["w1"]}, geometry=[Point(-106.5, 35.0)], crs="EPSG:4326") + gdf.to_parquet(root / "a" / "latest.parquet", write_covering_bbox=True) + res = _render(tmp_path, [_product("a")], source_root=str(root)) + assert res["a"]["providers"][0]["id_field"] == "id" + + +def test_bake_copies_the_chosen_file_and_serves_it_locally(tmp_path): + root, baked = tmp_path / "products", tmp_path / "baked" + _write_parquet(root, "pts", [_point("x:1")]) + _write_parquet(root, "poly", [_polygon("county")]) + (root / "poly" / "latest.geojson").write_text('{"type": "FeatureCollection", "features": []}') + res = _render( + tmp_path, [_product("pts"), _product("poly")], False, source_root=str(root), bake_to=baked + ) + assert res["pts"]["providers"][0]["data"]["source"] == f"{baked}/pts/latest.parquet" + assert res["poly"]["providers"][0]["data"]["source"] == f"{baked}/poly/latest.geojson" + assert sorted(str(p.relative_to(baked)) for p in baked.rglob("latest.*")) == [ + "poly/latest.geojson", + "pts/latest.parquet", + ] + + +def test_bake_fails_when_a_product_has_nothing_to_serve(tmp_path): + root = tmp_path / "products" + root.mkdir() + with pytest.raises(SystemExit, match="gone: no latest.geojson"): + _render(tmp_path, [_product("gone")], source_root=str(root), bake_to=tmp_path / "baked") diff --git a/tests/test_pygeoapi_patch.py b/tests/test_pygeoapi_patch.py new file mode 100644 index 0000000..a910621 --- /dev/null +++ b/tests/test_pygeoapi_patch.py @@ -0,0 +1,61 @@ +"""Tests for orchestration/pygeoapi/patch_pygeoapi.py.""" + +import importlib.util +from pathlib import Path + +import pytest + +PATCH = Path(__file__).resolve().parents[1] / "orchestration" / "pygeoapi" / "patch_pygeoapi.py" +spec = importlib.util.spec_from_file_location("patch_pygeoapi", PATCH) +pp = importlib.util.module_from_spec(spec) +spec.loader.exec_module(pp) + +# The lines patch_pygeoapi.py targets, as they appear in pygeoapi 0.25.dev0. +SOURCE = """ + if '/' in datetime_: + begin, end = datetime_.split('/') + if begin != '..': + begin = isoparse(begin) + if end != '..': + end = isoparse(end) + else: + target_time = isoparse(datetime_) + for batch in batches: + read += batch.num_rows + if read > limit: + batches_list.append(batch.slice(0, limit + 1)) + break + return gdf.__geo_interface__['features'][0] + number_matched = offset + len(rp) + result = gdf.__geo_interface__ +""" + + +def test_all_fixes_apply(): + out = pp.patch(SOURCE) + assert "limit + 1 - (read - batch.num_rows)" in out + assert out.count("_utc(isoparse(") == 3 + assert "def _utc(value):" in out + assert "__geo_interface__" not in out + assert out.count("show_bbox=False") == 2 + assert "number_matched = scanner.count_rows()" in out + + +def test_changed_source_fails(): + with pytest.raises(SystemExit, match="paging"): + pp.patch(SOURCE.replace("limit + 1))", "limit))")) + + +def test_already_patched_source_fails(): + with pytest.raises(SystemExit): + pp.patch(pp.patch(SOURCE)) + + +def test_utc_helper_treats_naive_as_utc(): + from datetime import datetime, timedelta, timezone + + ns = {} + exec(pp.UTC_HELPER, ns) + assert ns["_utc"](datetime(2020, 1, 1)).tzinfo == timezone.utc + mst = timezone(timedelta(hours=-7)) + assert ns["_utc"](datetime(2020, 1, 1, tzinfo=mst)).tzinfo == mst