Add configurable DEM sources and centerline file imports

.env.example

@@ -12,5 +12,30 @@ VITE_API_BASE=
 # Max number of stored REM job dirs before oldest are evicted (bounds disk use).
 MAX_COGS=200
 
+# Maximum GeoJSON, GeoPackage, or zipped-shapefile centerline upload (64 MiB).
+CENTERLINE_UPLOAD_MAX_BYTES=67108864
+# Basic processing limits protect public servers from unexpectedly complex files.
+CENTERLINE_MAX_FEATURES=10000
+CENTERLINE_MAX_VERTICES=500000
+
+# Optional custom DEM sources for the server engine. Both are disabled by
+# default, so public deployments do not gain upload/library access accidentally.
+DEM_UPLOAD_ENABLED=false
+# Maximum streamed upload size in bytes (default 20 GiB) and retention period.
+DEM_UPLOAD_MAX_BYTES=21474836480
+DEM_UPLOAD_TTL_HOURS=24
+# Expired uploads are removed at startup and on this repeating schedule.
+DEM_UPLOAD_CLEANUP_INTERVAL_MINUTES=15
+
+# To expose a read-only server DEM library, set the host directory and the
+# in-container path exactly as shown. Leave DEM_LIBRARY_DIR blank to hide it.
+DEM_LIBRARY_HOST_DIR=./dem-library
+DEM_LIBRARY_DIR=
+DEM_LIBRARY_LABEL=Server library
+
+# Maximum pixels retained from a custom DEM inside the current viewport.
+# Larger native-resolution windows are downsampled before RiverREM to bound RAM.
+CUSTOM_DEM_MAX_PIXELS=45000000
+
 # Short git commit hash shown in the footer (optional).
 GIT_SHA=dev

README.md

@@ -33,7 +33,7 @@ transparency) stays client-side and recolours instantly with no recompute.
 | | **Server** (default) | **Client** (beta) |
 |---|---|---|
 | Where | FastAPI + GDAL + RiverREM | entirely in the browser |
-| DEM | Mapterhorn tiles -> UTM mosaic | Mapterhorn tiles, read per-tile |
+| DEM | Mapterhorn, COG URL, streamed upload, or server library -> UTM | Mapterhorn tiles, read per-tile |
 | IDW | RiverREM: **k-NN, power 1**, KD-tree, `k` from sinuosity (faithful) | **all sampled points**, power configurable (QGIS-style) |
 | Output | single-band float32 COG (EPSG:3857) | `rem://` terrarium tiles built live, per tile |
 | Strength | fidelity, full-res, exports | no backend, live recompute, fast |
@@ -148,16 +148,69 @@ recolours instantly. The client engine skips all of that and runs in the browser
   "zoom": 13, "resolution_multiplier": 1, "idw_power": 2,
   "centerline_mode": "geojson",
   "centerline_geojson": { "type": "FeatureCollection", "features": [] },
-  "source_cog_url": null
+  "source_cog_url": null,
+  "source_dem_ref": null
 }
 ```
 -> `{ job_id, cog_url, dem_url, bounds:[w,s,e,n], rem_min, rem_max, width, height, river_name, river_length_m }`
 
-Other endpoints: `POST /cog/ingest` (reproject any-CRS COG to 3857), `POST /upload`
+Other endpoints: `GET /capabilities`, `POST /centerline/import`, `GET /dem/library`, `POST /dem/uploads`,
+`PUT /dem/uploads/{id}`, `POST /cog/ingest` (reproject any-CRS COG to 3857), `POST /upload`
 (zipped shapefile -> `upload_id`), `POST /thumb` (store a run thumbnail), `POST /runs/prune`.
 `source_cog_url` lets RiverREM use your own elevation COG (read via GDAL `/vsicurl/`)
 instead of Mapterhorn.
 
+Centerline uploads accept WGS84 GeoJSON, GeoPackage, and zipped shapefiles. GIS
+files with an assigned CRS are reprojected to WGS84 for the map; RiverREM then
+reprojects the normalized linework to the DEM CRS for computation. See
+[`docs/centerline-imports.md`](docs/centerline-imports.md).
+
+---
+
+## Using custom DEMs and centerlines
+
+Custom DEMs are available with the **Server** engine. Open **Compute**, choose a
+**DEM source**, and select one of:
+
+- **Mapterhorn** — downloads terrain for the current viewport; no setup required.
+- **COG URL** — enter an HTTP(S) URL for a remotely accessible elevation COG.
+- **Upload local DEM** — stream a georeferenced `.tif`/`.tiff` from the browser.
+  The server must enable uploads; uploaded DEMs expire and are removed automatically.
+- **Server library** — choose a DEM installed by the server owner. This option is
+  shown only when a library has been configured.
+
+To enable the optional sources, copy `.env.example` to `.env`, configure the
+desired settings, and rebuild the stack. A typical local configuration is:
+
+```dotenv
+DEM_UPLOAD_ENABLED=true
+DEM_UPLOAD_MAX_BYTES=21474836480
+DEM_UPLOAD_TTL_HOURS=24
+DEM_LIBRARY_HOST_DIR=/absolute/host/path/to/dems
+DEM_LIBRARY_DIR=/dem-library
+DEM_LIBRARY_LABEL=Available DEMs
+CUSTOM_DEM_MAX_PIXELS=45000000
+```
+
+The library is mounted read-only. On a public deployment it contains only DEMs
+installed by that server's administrator; leaving `DEM_LIBRARY_DIR` blank hides
+the library option. Large custom DEM viewports are reduced only when they exceed
+`CUSTOM_DEM_MAX_PIXELS`. The server checks that a custom DEM overlaps the current
+viewport and uploaded centerline before starting RiverREM. Full configuration,
+retention, and public-server guidance is in [`docs/dem-sources.md`](docs/dem-sources.md).
+
+To import a centerline, select **Draw**, click **Upload Centerline File**, and use:
+
+- `.geojson`/`.json` with EPSG:4326 longitude/latitude coordinates;
+- `.gpkg` with an assigned CRS (choose a layer if it contains several); or
+- `.zip` containing the `.shp`, `.shx`, and `.dbf` components of a Shapefile.
+
+GeoPackage and Shapefile linework is converted to EPSG:4326 for display. If its
+CRS metadata is missing, the UI asks for a source CRS such as `EPSG:6344`.
+The import summary shows the detected CRS, line/coordinate counts, length, and
+any skipped empty, invalid, or non-line features. See
+[`docs/centerline-imports.md`](docs/centerline-imports.md) for limits and validation details.
+
 ---
 
 ## Run it
@@ -182,7 +235,9 @@ docker compose -f docker-compose.local.yml up --build
 Open http://localhost:8088 (override with `WEB_PORT=NNNN`).
 
 Backend env: `TERRAIN_TILE_URL`, `TERRAIN_ENCODING` (`terrarium`|`mapbox`),
-`TERRAIN_MAX_ZOOM`, `PUBLIC_BASE`, `DATA_DIR`, `MAX_COGS`.
+`TERRAIN_MAX_ZOOM`, `PUBLIC_BASE`, `DATA_DIR`, `MAX_COGS`. Optional custom DEM
+uploads and a read-only server library are configured in
+[`docs/dem-sources.md`](docs/dem-sources.md).
 
 ---
 

backend/app/centerline.py

@@ -13,12 +13,175 @@
 from __future__ import annotations
 
 import os
+import json
+import math
 
 import geopandas as gpd
 import osmnx
+from osgeo import ogr
 from shapely.geometry import box, shape
 
 WATERWAY_TAGS = {"waterway": ["river", "stream", "tidal channel"]}
+CENTERLINE_MAX_FEATURES = max(1, int(os.environ.get("CENTERLINE_MAX_FEATURES", "10000")))
+CENTERLINE_MAX_VERTICES = max(2, int(os.environ.get("CENTERLINE_MAX_VERTICES", "500000")))
+
+
+class CenterlineCrsRequired(ValueError):
+    """Raised when a GIS dataset has linework but no assigned CRS."""
+
+
+def _dataset_layer_names(path: str) -> list[str]:
+    dataset = ogr.Open(path, 0)
+    if dataset is None:
+        raise ValueError("The centerline file could not be opened as a GIS dataset.")
+    names = [dataset.GetLayerByIndex(i).GetName() for i in range(dataset.GetLayerCount())]
+    dataset = None
+    return names
+
+
+def _crs_label(gdf: gpd.GeoDataFrame) -> str:
+    authority = gdf.crs.to_authority() if gdf.crs is not None else None
+    return f"{authority[0]}:{authority[1]}" if authority else str(gdf.crs)
+
+
+def _vertex_count(geometry) -> int:
+    if geometry.geom_type == "LineString":
+        return len(geometry.coords)
+    if geometry.geom_type == "MultiLineString":
+        return sum(len(line.coords) for line in geometry.geoms)
+    return 0
+
+
+def import_centerline_dataset(
+    path: str,
+    *,
+    display_name: str,
+    input_crs: str | None = None,
+    force_geojson_wgs84: bool = False,
+    layer_id_prefix: str = "",
+) -> list[dict]:
+    """Read all line layers and return WGS84 GeoJSON plus layer metadata.
+
+    RFC 7946 GeoJSON is always treated as EPSG:4326. Other formats use their
+    embedded CRS, or `input_crs` when the dataset does not declare one.
+    """
+    layer_names = _dataset_layer_names(path)
+    if not layer_names:
+        raise ValueError(f"{display_name} does not contain any GIS layers.")
+
+    results: list[dict] = []
+    missing_crs_layers: list[str] = []
+    unsupported_layers: list[str] = []
+    for layer_name in layer_names:
+        gdf = gpd.read_file(path, layer=layer_name)
+        if gdf.empty or "geometry" not in gdf:
+            continue
+        warnings: list[str] = []
+        usable_mask = gdf.geometry.notna() & ~gdf.geometry.is_empty
+        empty_count = int((~usable_mask).sum())
+        if empty_count:
+            warnings.append(f"Skipped {empty_count} empty geometr{'y' if empty_count == 1 else 'ies'}.")
+        gdf = gdf[usable_mask].copy()
+        line_mask = gdf.geometry.geom_type.isin(["LineString", "MultiLineString"])
+        unsupported_count = int((~line_mask).sum())
+        if unsupported_count:
+            warnings.append(
+                f"Skipped {unsupported_count} non-line feature{'s' if unsupported_count != 1 else ''}."
+            )
+        gdf = gdf[line_mask].copy()
+        if len(gdf) > CENTERLINE_MAX_FEATURES:
+            raise ValueError(
+                f"Layer '{layer_name}' contains {len(gdf):,} line features; "
+                f"the server limit is {CENTERLINE_MAX_FEATURES:,}."
+            )
+        # Inspect individual Shapely objects here rather than GeoSeries.length:
+        # source data can be geographic, where GeoPandas warns about interpreting
+        # length as a real-world distance (we only need a zero-length check).
+        nonzero_mask = gdf.geometry.map(lambda geometry: geometry.length > 0)
+        valid_mask = gdf.geometry.is_valid & nonzero_mask
+        invalid_count = int((~valid_mask).sum())
+        if invalid_count:
+            warnings.append(
+                f"Skipped {invalid_count} invalid or zero-length line{'s' if invalid_count != 1 else ''}."
+            )
+        gdf = gdf[valid_mask].copy()
+        if gdf.empty:
+            unsupported_layers.append(layer_name)
+            continue
+
+        vertex_count = sum(_vertex_count(geometry) for geometry in gdf.geometry)
+        if vertex_count > CENTERLINE_MAX_VERTICES:
+            raise ValueError(
+                f"Layer '{layer_name}' contains {vertex_count:,} coordinates; "
+                f"the server limit is {CENTERLINE_MAX_VERTICES:,}. Simplify the centerline and try again."
+            )
+
+        if force_geojson_wgs84:
+            # GeoJSON coordinates are defined as WGS84 longitude/latitude. Do
+            # not honor legacy/non-standard projected `crs` members silently.
+            gdf = gdf.set_crs("EPSG:4326", allow_override=True)
+            bounds = gdf.total_bounds
+            if (
+                not all(math.isfinite(float(value)) for value in bounds)
+                or bounds[0] < -180 or bounds[2] > 180
+                or bounds[1] < -90 or bounds[3] > 90
+            ):
+                raise ValueError(
+                    "GeoJSON coordinates must use WGS 84 (EPSG:4326) in "
+                    "[longitude, latitude] order."
+                )
+        elif gdf.crs is None:
+            if not input_crs:
+                missing_crs_layers.append(layer_name)
+                continue
+            try:
+                gdf = gdf.set_crs(input_crs, allow_override=True)
+            except Exception as exc:
+                raise ValueError(f"Invalid input CRS '{input_crs}': {exc}") from exc
+
+        source_crs = _crs_label(gdf)
+        try:
+            wgs84 = gdf.to_crs("EPSG:4326")
+        except Exception as exc:
+            raise ValueError(f"Could not convert layer '{layer_name}' to EPSG:4326: {exc}") from exc
+        bounds = wgs84.total_bounds
+        if (
+            not all(math.isfinite(float(value)) for value in bounds)
+            or bounds[0] < -180 or bounds[2] > 180
+            or bounds[1] < -90 or bounds[3] > 90
+        ):
+            raise ValueError(f"Layer '{layer_name}' produced coordinates outside the WGS84 range.")
+
+        try:
+            metric_crs = wgs84.estimate_utm_crs()
+            length_m = float(wgs84.to_crs(metric_crs).geometry.length.sum()) if metric_crs else 0.0
+        except Exception:
+            length_m = 0.0
+
+        layer_id = f"{layer_id_prefix}{layer_name}"
+        results.append({
+            "id": layer_id,
+            "name": layer_name,
+            "crs": source_crs,
+            "featureCount": int(len(wgs84)),
+            "vertexCount": vertex_count,
+            "lengthM": round(length_m, 1),
+            "warnings": warnings,
+            "geojson": json.loads(wgs84.to_json(drop_id=True)),
+        })
+
+    if missing_crs_layers:
+        raise CenterlineCrsRequired(
+            "CRS metadata is missing for line layer(s): " + ", ".join(missing_crs_layers)
+        )
+    if not results:
+        if unsupported_layers:
+            raise ValueError(
+                f"{display_name} has no usable centerlines. Centerlines must be non-empty "
+                "LineString or MultiLineString features."
+            )
+        raise ValueError(f"{display_name} does not contain a LineString or MultiLineString layer.")
+    return results
 
 
 def _features_from_bbox(west, south, east, north):
@@ -94,17 +257,25 @@ def geojson_to_shapefile(geojson: dict, out_path: str) -> str:
     # flatten to LineStrings, then merge touching ones
     lines = []
     for g in geoms:
-        if isinstance(g, LineString):
+        if isinstance(g, LineString) and not g.is_empty:
             lines.append(g)
         elif isinstance(g, MultiLineString):
-            lines.extend(g.geoms)
-    merged = linemerge(lines) if lines else None
+            lines.extend(line for line in g.geoms if not line.is_empty)
+    if not lines:
+        raise ValueError("The uploaded GeoJSON contains no usable LineString centerline.")
+    for line in lines:
+        for x, y, *_ in line.coords:
+            if abs(x) > 180 or abs(y) > 90:
+                raise ValueError(
+                    "Centerline coordinates must use WGS 84 (EPSG:4326) in longitude, latitude order."
+                )
+    merged = linemerge(lines)
     if isinstance(merged, LineString):
         out_geoms = [merged]
     elif isinstance(merged, MultiLineString):
         out_geoms = list(merged.geoms)
     else:
-        out_geoms = geoms
+        raise ValueError("The uploaded GeoJSON centerline could not be converted to line geometry.")
 
     gdf = gpd.GeoDataFrame({"name": ["centerline"] * len(out_geoms)}, geometry=out_geoms, crs="EPSG:4326")
     gdf.to_file(out_path)
@@ -113,7 +284,8 @@ def geojson_to_shapefile(geojson: dict, out_path: str) -> str:
 
 def normalize_uploaded_shapefile(upload_dir: str) -> str:
     """Find the .shp inside an uploaded/unzipped directory."""
-    for f in os.listdir(upload_dir):
-        if f.lower().endswith(".shp"):
-            return os.path.join(upload_dir, f)
+    for root, _dirs, files in os.walk(upload_dir):
+        for filename in files:
+            if filename.lower().endswith(".shp"):
+                return os.path.join(root, filename)
     raise FileNotFoundError("No .shp found in the uploaded shapefile bundle.")