|
| 1 | +# =============================================================================== |
| 2 | +# Copyright 2024 Jake Ross |
| 3 | +# |
| 4 | +# Licensed under the Apache License, Version 2.0 (the "License"); |
| 5 | +# you may not use this file except in compliance with the License. |
| 6 | +# You may obtain a copy of the License at |
| 7 | +# |
| 8 | +# http://www.apache.org/licenses/LICENSE-2.0 |
| 9 | +# =============================================================================== |
| 10 | +"""GeoPandas-backed persistence (see docs/framework-migration-plan.md, Phase A). |
| 11 | +
|
| 12 | +GeoPandas is the serialization engine for the product outputs. Every |
| 13 | +``dump_*_collection`` in ``ogc_features`` routes its features through a |
| 14 | +``GeoDataFrame`` — see :func:`route_feature_dicts_through_gdf`, called from |
| 15 | +``_dump_collection`` — so the OGC GeoJSON is emitted by GeoPandas rather than |
| 16 | +hand-assembled, and the same in-memory object can also write GeoPackage / |
| 17 | +PostGIS / GeoParquet. |
| 18 | +
|
| 19 | +The OGC FeatureCollection *envelope* (top-level id/title/timeStamp/links/ |
| 20 | +numberReturned and product-level extras) is a product concern that stays in |
| 21 | +``ogc_features._dump_collection``; this module only owns the per-feature geometry |
| 22 | ++ property serialization and the GeoParquet inter-asset handoff. |
| 23 | +""" |
| 24 | + |
| 25 | +import io |
| 26 | +import json |
| 27 | +from pathlib import Path |
| 28 | + |
| 29 | +import geopandas as gpd |
| 30 | +import pandas as pd |
| 31 | +from shapely.geometry import Point, shape |
| 32 | + |
| 33 | +from backend.constants import TDS |
| 34 | +from backend.persisters.ogc_features import _feature_id, _num, _tds_class |
| 35 | + |
| 36 | +# Property columns are the record's own keys minus the three that become |
| 37 | +# geometry. Kept here so the column order in the GeoDataFrame (and therefore the |
| 38 | +# GeoJSON `properties` order) matches the legacy feature output exactly. |
| 39 | +_GEOMETRY_KEYS = ("latitude", "longitude", "elevation") |
| 40 | + |
| 41 | + |
| 42 | +def _point(lat, lon, elev): |
| 43 | + """Mirror ogc_features._point_geometry: 2D Point unless an elevation is |
| 44 | + present, and no geometry at all when coordinates are missing (GeoPandas |
| 45 | + emits ``"geometry": null`` for those rows).""" |
| 46 | + if lat is None or lon is None: |
| 47 | + return None |
| 48 | + if elev is None: |
| 49 | + return Point(lon, lat) |
| 50 | + return Point(lon, lat, elev) |
| 51 | + |
| 52 | + |
| 53 | +def features_to_geodataframe(items) -> gpd.GeoDataFrame: |
| 54 | + """Build a GeoDataFrame from an iterable of ``(feature_id, geometry, props)`` |
| 55 | + where *geometry* is a shapely geometry or ``None`` and *props* is a dict. |
| 56 | +
|
| 57 | + The canonical persistence object. Property columns are held as ``object`` |
| 58 | + dtype so ``to_json`` preserves exact int/float/None types, and the feature id |
| 59 | + becomes the index (→ GeoJSON feature ``id`` via ``to_json(drop_id=False)``). |
| 60 | +
|
| 61 | + Note: a GeoDataFrame has a *uniform* column set, so a product whose features |
| 62 | + carry ragged property keys (e.g. per-well analyte pivots) gains explicit null |
| 63 | + columns here — a deliberate schema change. Uniform-schema products round-trip |
| 64 | + byte-identically. |
| 65 | + """ |
| 66 | + ids: list = [] |
| 67 | + geoms: list = [] |
| 68 | + rows: list[dict] = [] |
| 69 | + for feature_id, geom, props in items: |
| 70 | + ids.append(feature_id) |
| 71 | + geoms.append(geom) |
| 72 | + rows.append(props) |
| 73 | + frame = pd.DataFrame(rows, dtype=object) |
| 74 | + gdf = gpd.GeoDataFrame(frame, geometry=geoms, crs="EPSG:4326") |
| 75 | + gdf.index = ids |
| 76 | + return gdf |
| 77 | + |
| 78 | + |
| 79 | +def records_to_geodataframe(records: list) -> gpd.GeoDataFrame: |
| 80 | + """Build the canonical persistence GeoDataFrame from Summary/Site records. |
| 81 | +
|
| 82 | + - geometry: per-row Point (2D or 3D) in EPSG:4326, matching the legacy |
| 83 | + coordinate rule (elevation only when present). |
| 84 | + - columns: each record's ``keys`` except latitude/longitude/elevation, in |
| 85 | + order, plus ``tds_class`` for TDS records. |
| 86 | + - index: the OGC feature id ``"source:id"`` so ``to_json(drop_id=False)`` |
| 87 | + emits it as the feature-level ``id``. |
| 88 | +
|
| 89 | + Used for the GeoParquet inter-asset handoff (a source asset's summary/sites |
| 90 | + become a GeoDataFrame). Columns are held as ``object`` dtype so pandas does |
| 91 | + not coerce ints to floats or ``None`` to ``NaN`` — the emitted GeoJSON then |
| 92 | + preserves the exact value types the legacy hand-built features carried. |
| 93 | + """ |
| 94 | + rows: list[dict] = [] |
| 95 | + geoms: list = [] |
| 96 | + ids: list[str] = [] |
| 97 | + |
| 98 | + for r in records: |
| 99 | + props = {k: getattr(r, k) for k in r.keys if k not in _GEOMETRY_KEYS} |
| 100 | + if getattr(r, "parameter_name", None) == TDS: |
| 101 | + props["tds_class"] = _tds_class(_num(getattr(r, "latest_value", None))) |
| 102 | + rows.append(props) |
| 103 | + geoms.append( |
| 104 | + _point( |
| 105 | + getattr(r, "latitude", None), |
| 106 | + getattr(r, "longitude", None), |
| 107 | + getattr(r, "elevation", None), |
| 108 | + ) |
| 109 | + ) |
| 110 | + ids.append(_feature_id(getattr(r, "source", "") or "", getattr(r, "id", "") or "")) |
| 111 | + |
| 112 | + # object dtype preserves int/float/None exactly through to_json (see docstring). |
| 113 | + frame = pd.DataFrame(rows, dtype=object) |
| 114 | + gdf = gpd.GeoDataFrame(frame, geometry=geoms, crs="EPSG:4326") |
| 115 | + gdf.index = ids |
| 116 | + return gdf |
| 117 | + |
| 118 | + |
| 119 | +def geodataframe_to_features(gdf: gpd.GeoDataFrame) -> list[dict]: |
| 120 | + """Serialize a GeoDataFrame to a list of GeoJSON Feature dicts via GeoPandas. |
| 121 | +
|
| 122 | + ``drop_id=False`` promotes the index to each feature's top-level ``id``; |
| 123 | + ``na="null"`` emits missing values as JSON null (matching the legacy |
| 124 | + behavior of carrying ``None`` straight through).""" |
| 125 | + if gdf.empty: |
| 126 | + return [] |
| 127 | + return json.loads(gdf.to_json(drop_id=False, na="null"))["features"] |
| 128 | + |
| 129 | + |
| 130 | +def route_feature_dicts_through_gdf(features: list) -> list: |
| 131 | + """Rebuild a list of GeoJSON feature dicts by round-tripping them through a |
| 132 | + GeoDataFrame — the single hook that makes **every** ``dump_*`` in |
| 133 | + ``ogc_features`` GeoDataFrame-backed (called from ``_dump_collection``). |
| 134 | +
|
| 135 | + Geometry dicts are reconstructed with ``shapely.geometry.shape`` (point or |
| 136 | + polygon; ``None`` stays null). Because a GeoDataFrame has a uniform column |
| 137 | + set, products whose features carry ragged property keys gain explicit null |
| 138 | + columns (the chosen schema); uniform products round-trip byte-identically. |
| 139 | + Idempotent: features already emitted from a GeoDataFrame pass through |
| 140 | + unchanged. |
| 141 | + """ |
| 142 | + if not features: |
| 143 | + return features |
| 144 | + items = [] |
| 145 | + for f in features: |
| 146 | + geom_dict = f.get("geometry") |
| 147 | + geom = shape(geom_dict) if geom_dict else None |
| 148 | + items.append((f.get("id"), geom, f.get("properties", {}))) |
| 149 | + gdf = features_to_geodataframe(items) |
| 150 | + return geodataframe_to_features(gdf) |
| 151 | + |
| 152 | + |
| 153 | +def gdf_to_parquet_bytes(gdf: gpd.GeoDataFrame) -> bytes: |
| 154 | + """Serialize a GeoDataFrame to GeoParquet bytes for the Dagster inter-asset |
| 155 | + handoff (replaces the pickled ``_payload`` dicts). |
| 156 | +
|
| 157 | + The feature-id index is preserved (``index=True``) so it survives the |
| 158 | + round-trip and still becomes the OGC feature id downstream. GeoParquet keeps |
| 159 | + an explicit column schema + CRS in its metadata, so types round-trip far more |
| 160 | + faithfully than pickle-of-dicts — with one caveat: a column that mixes ints |
| 161 | + and ``None`` comes back float (Arrow encodes the null as a floating NaN). No |
| 162 | + summary/site column does that (``nrecords`` is always present), but the |
| 163 | + per-product fan-out should keep it in mind. |
| 164 | +
|
| 165 | + Requires the optional ``parquet`` extra (pyarrow); present in the |
| 166 | + orchestration deploy env via the dagster stack. |
| 167 | + """ |
| 168 | + buf = io.BytesIO() |
| 169 | + gdf.to_parquet(buf, index=True) |
| 170 | + return buf.getvalue() |
| 171 | + |
| 172 | + |
| 173 | +def parquet_bytes_to_gdf(data: bytes) -> gpd.GeoDataFrame: |
| 174 | + """Inverse of :func:`gdf_to_parquet_bytes` — read GeoParquet bytes back into a |
| 175 | + GeoDataFrame (geometry + CRS + feature-id index restored).""" |
| 176 | + return gpd.read_parquet(io.BytesIO(data)) |
| 177 | + |
| 178 | + |
| 179 | +# --------------------------------------------------------------------------- |
| 180 | +# Source-asset payload handoff (Parquet) — replaces the pickled {records, sites, |
| 181 | +# timeseries} dict crossing the Dagster IO manager. records/sites/observations |
| 182 | +# are flat scalar payload dicts (SummaryRecord/SiteRecord/ParameterRecord |
| 183 | +# _payloads), so a columnar Parquet table is the natural, typed handoff; the |
| 184 | +# combine rebuilds record objects from the dicts exactly as before. Geometry is |
| 185 | +# built later, in the dumpers — so these are plain Parquet, not GeoParquet. |
| 186 | +# --------------------------------------------------------------------------- |
| 187 | + |
| 188 | +# Marks which site-group each flattened observation belongs to, so the aligned |
| 189 | +# list-of-per-site-lists (sites[i] ↔ timeseries[i]) survives the round-trip. |
| 190 | +_SITE_IDX = "__site_idx" |
| 191 | + |
| 192 | + |
| 193 | +def _clean_nans(records: list[dict]) -> list[dict]: |
| 194 | + """pandas reads missing cells back as NaN; restore them to None so the |
| 195 | + rebuilt payload dicts match what was written (record classes expect None).""" |
| 196 | + return [{k: (None if pd.isna(v) else v) for k, v in row.items()} for row in records] |
| 197 | + |
| 198 | + |
| 199 | +def dicts_to_parquet_bytes(dicts: list[dict]) -> bytes: |
| 200 | + """Serialize a flat list of payload dicts (records or sites) to Parquet. |
| 201 | + object dtype keeps exact types; a uniform column set is fine — the record |
| 202 | + classes tolerate extra null-valued keys on rebuild.""" |
| 203 | + buf = io.BytesIO() |
| 204 | + pd.DataFrame(dicts, dtype=object).to_parquet(buf, index=False) |
| 205 | + return buf.getvalue() |
| 206 | + |
| 207 | + |
| 208 | +def parquet_bytes_to_dicts(data: bytes) -> list[dict]: |
| 209 | + """Inverse of :func:`dicts_to_parquet_bytes`.""" |
| 210 | + df = pd.read_parquet(io.BytesIO(data)) |
| 211 | + return _clean_nans(df.to_dict("records")) |
| 212 | + |
| 213 | + |
| 214 | +def timeseries_to_parquet_bytes(timeseries: list[list[dict]]) -> bytes: |
| 215 | + """Flatten the aligned list-of-per-site-lists to one Parquet table, tagging |
| 216 | + each observation with its site-group index so the grouping is recoverable.""" |
| 217 | + rows: list[dict] = [] |
| 218 | + for site_idx, site_ts in enumerate(timeseries): |
| 219 | + for obs in site_ts: |
| 220 | + rows.append({**obs, _SITE_IDX: site_idx}) |
| 221 | + buf = io.BytesIO() |
| 222 | + pd.DataFrame(rows, dtype=object).to_parquet(buf, index=False) |
| 223 | + return buf.getvalue() |
| 224 | + |
| 225 | + |
| 226 | +def parquet_bytes_to_timeseries(data: bytes) -> list[list[dict]]: |
| 227 | + """Inverse of :func:`timeseries_to_parquet_bytes` — regroup observations back |
| 228 | + into the aligned list-of-per-site-lists by ``__site_idx``. Group order is the |
| 229 | + site-group order (every persisted site carries ≥1 observation, so groups are |
| 230 | + contiguous 0..N-1).""" |
| 231 | + df = pd.read_parquet(io.BytesIO(data)) |
| 232 | + if df.empty: |
| 233 | + return [] |
| 234 | + n_groups = int(df[_SITE_IDX].max()) + 1 |
| 235 | + groups: list[list[dict]] = [[] for _ in range(n_groups)] |
| 236 | + for row in _clean_nans(df.to_dict("records")): |
| 237 | + idx = int(row.pop(_SITE_IDX)) |
| 238 | + groups[idx].append(row) |
| 239 | + return groups |
| 240 | + |
| 241 | + |
| 242 | +def write_geopackage(gdf: gpd.GeoDataFrame, path: str, layer: str) -> tuple: |
| 243 | + """Write a GeoDataFrame to a GeoPackage layer named *layer* and return its |
| 244 | + 2D bounds ``(minx, miny, maxx, maxy)`` in EPSG:4326. |
| 245 | +
|
| 246 | + The single GeoPackage writer for the GeoServer publish path. Sets a default |
| 247 | + CRS when absent, and flattens geometry to 2D — GeoServer's GeoPackage reader |
| 248 | + rejects a 3D CRS ("WGS 84 has 3 dimensions") when computing bounds, so |
| 249 | + elevation stays an attribute only. Bounds are returned so the caller can hand |
| 250 | + them to GeoServer explicitly (avoids a getBounds call that trips the same |
| 251 | + 3D-CRS bug). Raises on an empty frame (nothing to publish).""" |
| 252 | + if gdf.empty: |
| 253 | + raise ValueError(f"{layer}: GeoDataFrame has no features; nothing to write") |
| 254 | + if gdf.crs is None: |
| 255 | + gdf = gdf.set_crs("EPSG:4326") |
| 256 | + flat = gdf.copy() |
| 257 | + flat["geometry"] = flat.geometry.force_2d() |
| 258 | + minx, miny, maxx, maxy = (float(v) for v in flat.total_bounds) |
| 259 | + flat.to_file(path, driver="GPKG", layer=layer) |
| 260 | + return (minx, miny, maxx, maxy) |
| 261 | + |
| 262 | + |
| 263 | +def geojson_to_geopackage(geojson_path, layer_name: str, out_dir) -> tuple: |
| 264 | + """Convert a GeoJSON file to a GeoPackage whose layer is *layer_name* (so the |
| 265 | + published GeoServer layer takes that name). Returns ``(gpkg_path, bbox)`` with |
| 266 | + ``bbox`` the 2D EPSG:4326 bounds. Reads the GeoJSON with GeoPandas and writes |
| 267 | + the GPKG via :func:`write_geopackage`. Used by the GeoServer publish asset.""" |
| 268 | + gdf = gpd.read_file(geojson_path) |
| 269 | + gpkg_path = Path(out_dir) / f"{layer_name}.gpkg" |
| 270 | + bbox = write_geopackage(gdf, str(gpkg_path), layer_name) |
| 271 | + return gpkg_path, bbox |
0 commit comments