16b-1/2: terrain polygons as constraints (GeoPackage, GeoJSON, GML; CORINE maps) by skavhaug · Pull Request #108 · expertanalytics/rasputin

skavhaug · 2026-09-28T22:28:16Z

Increment 16b-1/2: terrain polygons (land cover, water, roads, rivers) as mesh constraints. Design: docs/increments/16b-terrain-polygons.md; acceptance: docs/benchmarks/2026-09-28/16b12-acceptance.md.

What it does

rasputin mesh --features PATH reads features from a GeoPackage (sqlite3 over its R-tree, no GDAL), GeoJSON, or the legacy CORINE GML (strict reader; the 2020 fixture repaired by tests/fixtures/corine/repair_gml.py).
CORINE class maps (corine, clc18_kode): borders become constraint edges carrying land_cover (bit 7) and, on the coast, water (bit 8). Shared borders are merged by the noder, with their bits combined.
Features are clipped to the domain as linework. A pre-clip in the source CRS drops only runs of edges that miss a 100 m margin (built in metres), never splits an edge, and widens the margin per edge for long edges (Ola, 2026-09-28).
features, features_crs, features_transform and features_notice (the Copernicus attribution) are recorded in .vtk/.ply. stderr reports features kept, dropped outside, clipped and empty skipped, plus a scan line for an unindexed GeoPackage.
Domain start chains moved from cli._domain_chains to chains.start_chains.

Numbers (battery, M1 Max)

A 287 km² catchment in EPSG:4326 over DTM10, at 1 m: 3.4 s without features, 6.3 s with the Norway extract, 14.0 s with the Europe-wide CORINE file. The two files give bit-identical meshes (1 235 226 triangles).
The quarter circle reproduces M5. tools/bench.py: no change against master (median +0.3 % over 42 cells; mesh sha256 identical).
48 km CORINE square at 1 m: 6.75 M triangles, node 0.15 s, refine 5.9 s.

Known open item

With the Europe file (EPSG:3035), features clip takes about 8.5 s: large far-away polygons are converted whole before cutting, because no bend bound is derived for 3035 to UTM yet. Follow-up agreed with Ola: derive the bound and cut before converting.

Size

Net +642 production lines (CLAUDE.md §2), under the 700 ceiling.

🤖 Generated with Claude Code

… features) Five suites (217 tests), support modules, the committed CORINE extract over the benchmark tile (tests/fixtures/corine/clc2018_7908_3.gpkg, 320 kB, EPSG:3035, cut by extract.py from Ola's GeoPackage) and a NOTICE covering it and the legacy GML (origin 8e30af4). test_features.py gains land_cover (7) and water (8); test_refine_golden.py takes the domain chains from chains.start_chains (R2's move). Pins in the increment file. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…ML, margin in metres, pre-clip keeps whole edges, I4 on the noded output, R10/R1/R4) Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…oints 2 and 3; 4 and 5 are corrections) Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…ixture, margin in metres, whole-edge pre-clip) test_io_gml.py: the reader is strict (GmlDocument has no `complete`; an unclosed root, a member left open, and both together are refused with GmlError naming the parser's line); the legacy file read is the repaired one, 200 Polygon features; tests of the repair script against the 2020 blob. tests/fixtures/corine/repair_gml.py adds the two missing closing tags (Ola: "Repear the broken file."), runs only on the 2020 file (git blob 97ddc98, 8e30af4), and repaired the committed file; NOTICE entry 2 says so. test_feature_input.py: source_region is the domain plus about 100 m in the DEM's CRS, for 4326, 3035 and 25833 sources; pre_clip keeps whole edges (closed ring kept with its start, one run joined across the start, two runs, a line losing its middle or ends, a corner touch, no vertex added, nothing kept); a 10 km EPSG:4326 east-west edge across the region boundary gives the no-pre-clip noded edges (an intersecting pre-clip moves the domain crossing 3.26 m). The increment file pins the two new names. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…ures

…16d (a record of the input polylines, proposed by Ola) added; 16b-0 shipped with #107 Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…imits per persona (Ola) Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

… dropped (Ola) Ola's two rulings of 2026-09-28, recorded in "Ruled by Ola", R5, R10, I4, the LOC table and the red-suite consequences: - "yes, widen per edge": each edge's pre-clip margin is widened by a bound on how far its source-straight and DEM-straight lines separate. A stated bound for a geographic source over a Transverse Mercator DEM (checked with pyproj), w = 0 for the DEM's own CRS, and for any other pair (EPSG:3035 included) the feature is reprojected first and pre-clipped exactly in the DEM's CRS. The GeoPackage query is widened per row, and the GeoJSON/GML intersects filter is removed. I4 holds with no length limit. - "A table in a run report is output, not data.": no edge_vocabulary field and no report line; the .vtk and .ply vocabulary records already exist. Also fixes the stale Status line (Q1-Q6 are ruled) and the GML feature count in "Found in passing" (200, not 399). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

For design 1e67a72 (R5 "Long edges", R10): - test_cli_mesh_features.py: test_the_edge_vocabulary replaced by a test on the writer's feature_bits/feature_names (7 land_cover, 8 water) and feature_vocabulary, with no edge_vocabulary field; test_the_edge_vocabulary_without_features and the docstring item removed. - test_feature_input.py: a 60 km EPSG:4326 edge at 70 N over UTM 33, ~40 m outside the region, whose UTM-straight line crosses the domain ~194 m north of its midpoint, is kept and gives the no-pre-clip noded edges (I4); the test shows first that it misses the region. One unit test pins pre_clip(geometry, region, widening) with widening(a, b) per edge and a closed distance test. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…le-edge pre-clip, CLI --features Design: docs/increments/16b-terrain-polygons.md (rulings d84701f, 413e91b, 1e67a72). Red suite: e99c8ea, 3990449, 972312c. No test file touched. - io/geopackage.py: header decode (every envelope code, either byte order, empty and extended flags), layer_info, query_features in primary-key order; with an R-tree, each row's box widened by s*max(1, s/H) (R5, "Long edges"). - io/gml.py: strict GML2 reader over ElementTree; refusals name the line. - io/repository.py: open_geopackage, read-only URI; still the only io module that opens files. - feature_input.py: class maps (property, corine, corine-water, clc18_kode), sources by suffix, region in metres in the DEM's CRS, pre_clip(geometry, region, widening) keeping whole edges, w(e) for a geographic source over a Transverse Mercator DEM, reproject-first with w = 0 for any other pair, line clip to the domain. - chains.py: start_chains; cli._domain_chains moved here. - features.py: land_cover (bit 7) and water (bit 8). - crs.py: transform_definition, to see grid steps in PROJ's pipeline. - cli.py: --features, --features-crs, --features-layer, --features-map; fields features, features_crs, features_transform, features_notice; stderr counts; --stats rows "features read" and "features clip". Three tests stay red, each a test defect reported to the orchestrator: test_a_ring_crossing_the_boundary_is_open_pieces (asserts 200 == 250), test_stderr_and_stats (^-anchored regex over whitespace-collapsed output), test_the_uniform_mesh_is_unchanged (file digest predates Q3's two bits). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

- test_a_ring_crossing_the_boundary_is_open_pieces: the literal 250 was wrong; the ring's length inside x <= 300 is 50 + 100 + 50 = 200, which the test's own shapely oracle on the line above already says. - test_stderr_and_stats: invoke() collapses output onto one line, so the MULTILINE ^ anchor could never match; the rows are found by their cell, "| features read |". Checked to fail with the row removed. - test_the_uniform_mesh_is_unchanged: re-recorded. Compared with 972312c's output, only feature_bits, feature_names and feature_vocabulary differ (bits 7 land_cover and 8 water added); points, cells and all other fields are the same. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>