diff --git a/.github/workflows/main.yaml b/.github/workflows/main.yaml index 26e89e62..3cc62803 100644 --- a/.github/workflows/main.yaml +++ b/.github/workflows/main.yaml @@ -185,15 +185,20 @@ jobs: - name: Tests run: pytest - # The GeoTIFF suites branch on the `codecs` extra: some tests need - # imagecodecs (LZW, the real Kartverket DEM), others test its absence. - # No single environment runs both halves, so the step above covers the - # absent half and this one installs the extra and re-runs both files. + # Suites that branch on the `codecs` extra: some tests need imagecodecs + # (LZW, the real Kartverket DEM and DTM10 extracts), others test its + # absence. No single environment runs both halves, so the step above + # covers the absent half and this one installs the extra and re-runs the + # codecs-gated suites: the GeoTIFF ones, and the CLI and golden suites that + # mesh the real tile (test_cli_mesh_dem.py runs in the viewer step below). - name: Tests with the codecs extra run: | python -m pip install -e ".[dev,codecs]" pytest --no-cov tests/python/test_io_geotiff.py tests/python/test_geotiff_fixtures.py \ - tests/python/test_io_read_meta.py tests/python/test_cli_mesh_mosaic.py + tests/python/test_io_read_meta.py tests/python/test_cli_mesh_mosaic.py \ + tests/python/test_cli_mesh_domain_crs.py tests/python/test_cli_mesh_domain.py \ + tests/python/test_refine_golden.py tests/python/test_cli_constraint_feet.py \ + tests/python/test_cli_mesh_refine.py tests/python/test_cli_start_quality.py # Increment 13, ruling 10: read the .vtk output back with the readers # ParaView uses. The suite skips when vtk is absent, so it is run only diff --git a/ROADMAP.md b/ROADMAP.md index 09dc9062..e5eed1f1 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -44,7 +44,7 @@ increment that most needs a picture to check against | — | `tools/bench.py`: the 1 m benchmark and the thread-scaling sweep from one checked-in command, with power state, quality and commit recorded per run; the one-off scripts in `docs/benchmarks/2026-09-26/` are its specification. Rule 2's acceptance run needs it | shipped with branch `tools-bench`'s PR | `docs/benchmarks/bench-py.md` | | — | The serial phase: profile refine's serial insert-and-flip phase, then parallelise what the profile blames. Scaling tops out at about 2.0-2.2×. Profiled 2026-09-27: serial part about a third of 1-thread refine, mostly Lawson legalisation; the scan stops near 5× from load imbalance | profiled; designed as increment 21 (Ola's rulings 2026-09-27: L1 determinism, one path, at most 2 % more triangles). **21a shipped with branch `increment21a-quick-wins`'s PR**: dynamic scan blocks, active merge, reused flip stack; mesh bit-identical; on AC refine -10 to -15 % at 8 threads, ceiling 2.1x -> 2.3x (`docs/benchmarks/2026-09-27/21a-acceptance.md`). **21b shipped with branch `increment21b-lattice-incircle`'s PR**: an int64 lattice incircle answers 99.97-100 % of refine's incircle tests; mesh bit-identical; on battery refine -15 % at 8 threads, ceiling 2.3x -> 2.5x (`docs/benchmarks/2026-09-27/21b-acceptance.md`). **21c measured** (`docs/benchmarks/2026-09-27/21c/README.md`): option C costs 6.5-10 % more triangles; A1 about 1 %; A0 is bit-identical (on these two inputs, not proven), and with evaluate-once is modelled at only 7-11 % faster refine at 8 threads and slower at 4. **21d is deferred** (Ola, 2026-09-27: "Review and push 21c, then basin-work"), behind the work in "Order of work" below; domain decomposition parallelises scan and split together | `docs/increments/21-parallel-refine.md`, `docs/benchmarks/2026-09-27/serial-profile/README.md` | | — | Release hardening: measure libc++'s `_LIBCPP_HARDENING_MODE_FAST` (and libstdc++'s assertions on the GCC leg), which bounds-check `std::vector`, `span` and the like in Release, on the 1 m benchmark; switch it on if the cost is small. Today an out-of-range index in shipped code the tests miss is undefined behaviour and crashes Python (Ola, 2026-09-27: "I'm surprised we don't have proper memory control"; CI's ASan/UBSan/TSan cover what the tests reach) | to measure (Ola, 2026-09-27); 21d, which it was placed after, is deferred behind the basin work | none yet | -| 15 | A DEM in several tiles, and the domain in its own CRS: 15a Norway, many tiles in one projected CRS (Ola's 254 DTM10 UTM33 tiles; only the selected tiles must share a lattice, since 8 of them sit half a cell off); 15b the domain polygon reprojected into the DEM's CRS; 15c and 15d the São Francisco basin (geographic DEMs meshed in the DEM's own lattice frame, not resampled; window decoding and memory) | **15a shipped with branch `increment15-dem-mosaic`'s PR**: `--dem DIR` or several files, `--bbox`, tiles selected and stitched on one lattice, Q2-Q5 refusals; overlaps that disagree (real DTM10 tiles exported on different dates, up to 52 m) are split down the middle and reported per seam (Ola's Q1 revised, 2026-09-28). The design's acceptance box, 9 tiles and 10,051² nodes, meshes: 11.05 M triangles in 31 s on battery. 15b next (Ola, 2026-09-27: Norway first) | `docs/increments/15-dem-mosaic.md` | +| 15 | A DEM in several tiles, and the domain in its own CRS: 15a Norway, many tiles in one projected CRS (Ola's 254 DTM10 UTM33 tiles; only the selected tiles must share a lattice, since 8 of them sit half a cell off); 15b the domain polygon reprojected into the DEM's CRS; 15c and 15d the São Francisco basin (geographic DEMs meshed in the DEM's own lattice frame, not resampled; window decoding and memory) | **15a shipped with branch `increment15-dem-mosaic`'s PR**: `--dem DIR` or several files, `--bbox`, tiles selected and stitched on one lattice, Q2-Q5 refusals; overlaps that disagree (real DTM10 tiles exported on different dates, up to 52 m) are split down the middle and reported per seam (Ola's Q1 revised, 2026-09-28). The design's acceptance box, 9 tiles and 10,051² nodes, meshes: 11.05 M triangles in 31 s on battery. **15b built with branch `increment15b-domain-crs`'s PR**: `--domain` in its own CRS (e.g. EPSG:4326) over a UTM33 mosaic, reprojected vertex by vertex. Then 16b (Ola, 2026-09-27: Norway first) | `docs/increments/15-dem-mosaic.md` | | 16b | Interior polygons and polylines as constraints ("terrain polygons": lakes, land cover, roads, rivers): `--features PATH`, a GeoJSON `FeatureCollection`, each feature naming a vocabulary property; closed or open `Breakline`s, crossings noded, off-node vertices with bilinear z. **Its working example is real data** (Ola, 2026-09-27): CORINE Land Cover 2018 over the benchmark tile `7908_3_10m_z33.tif`. The source is Ola's local copy, `rasputin_data/corine_sql/.../U2018_CLC2018_V2020_20u1.gpkg` (8.2 GB, EPSG:3035, a sibling of this repository, which also holds the 254-tile DTM10 archive for gap 6). It reads without GDAL: sqlite3 over its R-tree, the GeoPackage blob header stripped, `shapely.wkb`, then `pyproj` to 25833, so CRS stops in Python as before. Probed 2026-09-27 in 0.2 s: 60 polygons in 8 classes (heath, bare rock, sparse vegetation, bogs, intertidal flats, water, sea, urban), 11 068 vertices clipped to the tile, median segment 54 m against 10 m cells. The EEA's public ArcGIS service (`image.discomap.eea.europa.eu`, `Corine/CLC2018_WM`) returns the same 11 068 clipped vertices and is the route for anyone without the file. What it forces on 16b's design: neighbouring polygons share their boundaries, so each shared edge arrives twice; the polygons run past the domain and must be clipped; the extract is committed as a fixture with the Copernicus attribution. Placed before 20c because 20c may split constraint segments and should be designed and measured on inputs that have interior ones | designed in 16's R6, no increment file yet; after the serial phase, before 20c | `docs/increments/16-domain-polygon.md` (R6) | | 20c | Soft quality criterion: a penalty that each Steiner node or constraint split must pay for in angle gained, instead of 20's hard 25°; applied at the start and during DEM refinement; may split constraint segments when that improves the mesh. Ola's rulings on 20's C1-C3 | to design after 16b (`@architect` measures cost against 20 first) | `docs/increments/20-start-quality.md` (Ola's rulings) | | — | Auto-catchment: the watershed upstream of a coordinate, computed from the DEM and handed to `--domain`, so a catchment no longer has to be supplied as a file (Ola, 2026-09-27: "not far into the future"). The textbook route is depression handling (Priority-Flood, Barnes, Lehman and Mulla 2014), D8 flow directions (O'Callaghan and Mark 1984) and accumulation, the pour point snapped to the strongest flow nearby, the upstream cells traced and their outline turned into a polygon; the literature check is `@architect`'s. Open for its design: whether it runs in the C++ core (a 10 m tile is 25 M cells); how a stair-stepped cell outline becomes a domain polygon, which meets input coarsening; and that a real catchment crosses tile edges, so it needs gap 6 (a DEM in several tiles) first. Legacy has nothing on it (`grep -rliE "watershed|flow.?acc|flow.?dir|pour.?point|catchment" legacy` returns no files) | to design; after gap 6, which it needs; placed after 20c, can move ahead of it on Ola's word | none yet | @@ -100,9 +100,10 @@ Norwegian cases sorted first"): 15a and 15b (Norwegian multi-tile DEM), then Open and not MVP-blocking: **inputs in their own CRS**. The user (2026-09-26): "I don't think the domain CRS should have to match the DEM CRS in the future. This must be written down. The questions will be in which coordinate system we -shall do the math." Increment 16 requires a match for now; a later increment -chooses the computation CRS and reprojects inputs into it (see -`docs/increments/16-domain-polygon.md`, "Ruled by the user"). Also **a general point insertion policy**, which the user +shall do the math." **Settled for the domain and a projected DEM by increment +15b**: the domain comes in its own CRS and is reprojected into the DEM's, which +is the computation CRS. Feature geometry (16b) follows the same path; a +geographic DEM's computation frame is 15c (`docs/increments/15-dem-mosaic.md`). Also **a general point insertion policy**, which the user wants to discuss (2026-09-26): which points refinement inserts, beyond increment 14b's worst DEM node, and how that meets the sizing field, coastline constraints and points that are not DEM nodes. The user's direction for it @@ -115,6 +116,16 @@ coarsening input geometry (`parallel_refinement.md` step 2), 5d above, and land-cover partitioning, whose foundation is increment 7's property sets and whose consumer does not exist. +Open, for later (Ola, 2026-09-28): **the surface model, DOM10, alongside the +terrain model.** Kartverket's DOM10 is the top surface (tree crowns, roofs, +bridges) where DTM10 is the bare ground; both are published on hoydedata.no +and, where laser coverage exists, derive from the same NDH laser data (DTM10 +falls back to the 2013 contour model elsewhere). Two uses Ola wants to keep open: **shading** (terrain and surface +shadowing, e.g. for solar radiation), which needs the surface, not the ground; +and **canopy and building height, DOM10 minus DTM10**, for vegetation-related +work and land-cover classification. Hydrology keeps meshing the DTM: a DOM +mesh would dam rivers at bridges. Not designed; no increment yet. + Known defect: increment 14's NoData carving from a NoData corner appears to advance one node per round. Meshing the real tile from its outline alone took 5 054 rounds and 109 s against 0.35 s from the stride grid diff --git a/docs/increments/15-dem-mosaic.md b/docs/increments/15-dem-mosaic.md index 8736b157..209298d8 100644 --- a/docs/increments/15-dem-mosaic.md +++ b/docs/increments/15-dem-mosaic.md @@ -1,8 +1,9 @@ # Increment 15 — a DEM in many tiles, inputs in their own CRS, and the computation frame Status: **Q1-Q5 ruled by Ola; 15a implemented** on branch -`increment15-dem-mosaic` (red `2696bc2`, green `ff7cc8d`), in review. 15b-15d -are designed, not implemented; Q6-Q10 are open. Written by `@architect` before +`increment15-dem-mosaic` (red `2696bc2`, green `ff7cc8d`; merged as `40de334`, #105). **15b +implemented** on `increment15b-domain-crs` (red `372bb99`, green `0011741`). +15c-15d are designed, not implemented; Q6-Q10 are open. Written by `@architect` before `@tester`, per `docs/increments/README.md` step 1. ## Ruled by Ola @@ -58,12 +59,28 @@ are designed, not implemented; Q6-Q10 are open. Written by `@architect` before Ola will also download a fresh single-date set from hoydedata.no, whose overlaps should mostly agree (the probe found one same-date pair that does not); the fresh set is not blocking ("the data is data", Ola). + **Measured on arrival (2026-09-28, `DTM10_UTM33_20260925`):** the same 254 + tiles, format and lattices (the eight half-cell tiles remain), and still + exported tile by tile (side-file dates 2020-06 to 2026-09); 144 tiles were + re-exported since the 2022 download, and the other 110 are byte-identical. + The probe (120 sampled pairs, 119 overlapping, seed 1): 32 of 32 same-date pairs agree within 1 mm; + 34 of 87 different-date pairs do not, worst 34.1 m (7404_2 | 7404_3). The + 15a acceptance box reports 5 seams (4 before), up to 22.3 m. So the fresh + set does not remove disagreeing overlaps; the midline rule and the seam + report stay necessary. - **2026-09-28, the seam report ignores differences below 1 mm** (Ola: "Ignore below 1mm"). A seam counts, and the report lists, only nodes where |a − b| ≥ 1 mm; float noise such as 7807_1 | 7808_4 (31 nodes, 1.5e-5 m) and 7910_2 | 7910_3 (4 nodes, 2.4e-7 m) no longer appears. The midline rule itself is unchanged: which tile's value a node takes does not depend on the threshold. +- **2026-09-28, with `--domain` the seam report counts only nodes inside the + needed region** (Ola: "yes, go with a"). The needed region is the domain + grown by the chosen plan's cell diagonal (15b), exactly the nodes bilinear z + can read. A seam whose overlap lies inside the plan's rectangle but wholly + outside the needed region is no longer reported (an irregular catchment's + bounding box can be several times its area). Without `--domain` the report + is unchanged. Which tile's value a node takes is unchanged everywhere. - **2026-09-27, Q5 as read after review ("Yes", to the main session's proposal):** @reviewer found that DTM10's 51-node overlaps make a box that only reaches into a half-cell tile's overlap strip select that tile and be @@ -679,11 +696,11 @@ before anything else sees them. straight segment between them in the frame the mesh is built in. Within one UTM zone the edges bend by millimetres per kilometre of edge (B7's order of magnitude). Not densified. -- `check_crs`'s must-match rule (`domain.py:97`) is replaced by the transform. - It was kept as "one replaceable function at the Python boundary" for exactly - this (16, "Ruled by the user"). The extent check (16 R1) now runs after the - transform, against the mosaic's coverage (R4 point 5), not one tile's - rectangle. +- `check_crs`'s must-match rule (`domain.py:97` at `d34d79d`) is replaced by + the transform. It was kept as "one replaceable function at the Python + boundary" for exactly this (16, "Ruled by the user"). The extent check (16 + R1) now runs after the transform, against the mosaic's coverage (R4 point + 5), not one tile's rectangle. - A GeoJSON file without a `crs` member is EPSG:4326 (RFC 7946), and is now transformed rather than refused. - **Provenance:** the `.vtk` records `domain_crs` (the input's CRS) and @@ -858,6 +875,9 @@ order. 15c and 15d follow later. | | `dem_input.py`: `DemRequest`, `open_dem` | 35 | | | | `cli.py`: `--dem` list, `--bbox`, refusals, fields | 45 | | | | **15a total** | **375** | **520** | +| | *15b as built, measured at review (d34d79d..914dfc8, CLAUDE.md §2 unit): 150 added, of which `dem_input.py` 48 against 20 estimated; 15a + 15b against master 675, under the ceiling; they ship as two PRs per Ola's order* | | | +| | *15b after the merge of 15a, against 15a's tip (`201a4e7..80fda64`): 148 added, 97 net* | | | +| | *15b at `e6b69de`, against 15a's tip (`201a4e7..e6b69de`): 162 added, 108 net; 2 over the 160 worst case, the overrun being the unestimated seam mask for Ola's 2026-09-28 ruling (`mosaic.py` +13); well under the ceiling, no split* | | | | | *15a as built, measured at review (branch diff, CLAUDE.md §2 unit): 533, of which `mosaic.py` 325 against 200 estimated* | | | | | *15a after Ola's Q1 revised (cba0073): 596 added against master, of which `mosaic.py` 372 against 200 estimated* | | | | | *15a at the tip (a253a77): 597 added against master, blank lines excluded (`CLAUDE.md` §2 as ruled by Ola 2026-09-28); 554 net* | | | @@ -1196,6 +1216,155 @@ for dataset `dddbb667-1303-4ac5-8640-7ec04c0e3918`: "Åpne data", CC BY 4.0): Deflate, 0.6 MB together. +### Pinned by the red suite (15b) + +Names and behaviours the design left open, fixed by the red commit's suites: +`tests/python/test_crs.py`, `test_domain.py` (rewritten for the 15b API), +`test_dem_input_domain.py` and `test_cli_mesh_domain_crs.py`. 15b is an +ordinary suite (no invariant-critical suite, so no full mutation round), but +the axis-order mutants were run: `test_crs.py` against a scratch `crs.py` +(not committed) passed 31 of 31, and killed `always_xy=False`, `always_xy` +dropped, input columns swapped, output columns swapped, pyproj's `CRSError` +let through, and a float32 round trip. Nothing else was run against a +scratch implementation. + +**`tin_engine/crs.py`.** + +- `parse_crs(text) -> pyproj.CRS`: anything `CRS.from_user_input` accepts, + including WKT2 and a PROJ string with no EPSG code. Anything else is a + `ValueError` naming the text. pyproj's `CRSError` is a `RuntimeError`, so it + must not leak: every `except ValueError` that makes a usage error would miss it. +- `reprojector(src, dst)` takes text or `pyproj.CRS` for each, and returns a + callable from an `(N, 2)` array-like of `(x, y)` to an `(N, 2)` float64 + array, bit for bit pyproj's `always_xy` transform. `x` is easting or + longitude whatever the CRS's axis order. +- The only `from_crs` call in `src_python/` is in `tin_engine/crs.py` (an AST + scan, whose finder is tested on planted source). + +**`tin_engine/domain.py`.** + +- `read_domain(path, crs=None) -> DomainPolygon`. It no longer takes the DEM: + the domain is read before the tiles are planned, since its bounds choose + them. Parsing, the geometry refusals and the file-against-`--domain-crs` + refusal are 16's, unchanged; a CRS with no EPSG code and a GeoJSON `crs` + member naming OGC's CRS84 are now read. `check_crs` is gone. +- `DomainPolygon`: `polygon` and `crs: str` (text pyproj parses to the + domain's CRS); no `epsg` field. +- `DomainPolygon.to_crs(dst)` (text or `pyproj.CRS`): + - every vertex of every ring is pyproj's `always_xy` transform of the input + vertex, bit for bit; the vertex count per ring is unchanged; the source is + unchanged; + - the result is re-oriented to the winding contract (outer counter-clockwise, + holes clockwise). A CRS whose easting points west (`+axis=wnu`) mirrors the + ring, and is the test; + - **same CRS is decided by `pyproj.CRS` equality, not by `to_epsg()`**. UTM 33 + with a west-pointing axis answers `to_epsg() == 25833` at pyproj's default + confidence, and skipping its transform would put the domain 1000 km west. + For an equal CRS the coordinates are returned bit for bit and no + transformer is made (`Transformer.from_crs` is patched to fail); + - a vertex with no image in `dst` (pyproj answers `inf`, as for latitude 95 + or a UTM easting read as a longitude) is a `DomainError` naming both CRSs, + raised at reading or at the transform. +- `check_extent(domain, meta)`: 16's node-rectangle check as a public + function, closed on the border, holes included; the message says `outside` + and names the vertex. + +**`tin_engine/dem_input.py`.** + +- `DemRequest(sources=, bounds=, nodata=, domain=None)`, `domain` a + `DomainPolygon` as read. Both `bounds` and `domain` is a `ValueError`. +- `open_dem` moves the domain into the DEM's CRS, plans on its bounds with its + needed region, runs `check_extent` against the plan's `meta`, and only then + assembles. `DemInput.domain` is the domain in the DEM's CRS, `None` without + one. The plan equals the plan of `bounds` equal to the moved domain's bounds, + except where a domain vertex lies within the 1e-6-cell snap band past a node + line: there `_past` moves that edge out and the domain's window is one node + line wider than `--bbox`'s (found at green, 15b). On one file the plan is a + window of it. +- **The snap band** (review S1, `TestTheSnapBand`): with every edge on a node + line, or one edge 2e-6 cell past one, the domain's plan equals `--bbox`'s; + with one edge 1e-7 cell past a node line it is one node line wider on that + side only, for each of the four edges. A vertex exactly on a node line + widens nothing. +- **The needed region**, "the domain polygon grown by one cell", is pinned + only away from its edge: a missing node 0.73 cell from a domain vertex is + refused (`in no tile`, naming its x); missing nodes 2.5 cells or more from + the domain, including under a hole of the domain, are NaN filler; a domain + enclosing a missing tile is refused (the interior is needed). A growth of + exactly one cell, and the metric (Euclidean or per axis), are not ruled. +- **The cell is the chosen plan's** (review B1, + `TestNeededRegionIsGrownByThePlansSpacing`): in a repository holding two + spacings in one EPSG, a tile the domain does not select, and whether its + name sorts first or last, changes neither the plan nor the refusal. A + missing node inside one 10 m cell but outside one 1 m cell is refused on a + 10 m plan and NaN filler on a 1 m plan. When growing the region changes the + chosen lattice, the region is grown again by the new lattice's cell (the + re-plan loop; @reviewer showed a once-only growth passed this class). +- **Open, recorded at review (not ruled):** (a) when the ungrown domain already + misses a node, the refusal counts only the domain's own missing nodes, a + lower bound on the grown region's; (b) with several spacings in one EPSG + (at least three lattices, @reviewer's reasoning, not a run) the lattice + choice can alternate as the region grows; the loop settles it by keeping the + larger growth, which can only add refusals. Not reachable on Ola's DTM10 + archive (one spacing). (c) `_past`'s snap-band re-plan does not re-run the + growth loop; if its half-cell widening selected a lattice with a bigger + cell, the smaller growth would be kept. Only inside the 1e-6-cell band with + several lattices. +- **One CRS** (review S2, `TestOneCrs`): tiles in more than one EPSG code with + a domain are a `MosaicError` before any `load`, naming the count, the codes + and "a domain needs one". +- Every extent refusal (no tile, uncovered, outside) fires before any `load`. +- **Seams inside the needed region only** (Ola's ruling of 2026-09-28; test + amendment, `TestSeamsInsideTheNeededRegion`, and + `test_cli_mesh_domain_crs.py::TestSeamsWithADomain::test_only_the_needed_region_is_counted`). + On 21 x 21 quadrant tiles (dx = dy = 10) with `ne.tif` planted off `nw.tif` + at five nodes of their shared column: a thin diagonal strip whose plan is + the whole grid, and whose region misses every planted node, reports + `seams == ()` and `dem_seams` `none`; an L whose region takes two of them + reports exactly those (nodes 2, max 2, median 1.25), equal to `seams_of` + masked by the domain grown by the plan's cell diagonal, mitred + (`mosaic_fixtures.seams_of` takes an optional node mask); a rectangle 7 m + short of the column counts the nodes 7 m outside it. `--bbox` at each + domain's bounds has the same plan and reports all four planted nodes, and + the domain's mosaic is `--bbox`'s bit for bit and the midline oracle's. + Red at `5d12ad0` for the strip and the L (both files); the other cases pass + there as pins. Against a scratch implementation (a node mask in `assemble`'s + seam loop, not committed): masking by the ungrown polygon fails the L and + the 7 m case, no masking fails the strip and the L. + +**`cli.py`.** + +- `--bbox` with `--domain` is a usage error naming both. +- `domain_crs` is `EPSG:n` when pyproj finds an **exact** EPSG code (so OGC's + CRS84 is not recorded as `EPSG:4326`), and otherwise ASCII text pyproj parses + back to an equal CRS. `domain_transform` is the `description` of + `Transformer.from_crs(domain CRS, DEM CRS, always_xy=True)`, and contains + `none` (any case) for a domain already in the DEM's CRS. Only the `.vtk` + fields are pinned, not a `.ply` comment. +- A transformed domain refused for its extent (no image, no tile, outside) is + a usage error naming the domain's CRS and the DEM's EPSG code, and writes + nothing; one reaching past the tiles also says `--domain` and `outside`. +- **Same CRS, bit for bit.** `test_cli_mesh_domain_crs.py`, relational on one + machine: for the micro-TIFF square of 16's suite at 1 m and the quarter + circle on the benchmark tile at 10 m (`needs_codecs`), the run equals (SHA-256 + of points, cells, cell and point arrays; not field data) the same run with + `cli.open_dem` replaced by 16's data flow, the whole file and the domain as + `read_domain` returned it; the domain `_dem_mesh` receives has the read + vertices bit for bit; `Transformer.from_crs` is never called. It replaced + digests recorded at `d34d79d` on macOS arm64, whose square Linux x86 (GCC) + does not reproduce (PR #106's CI). The platform-stable absolute anchor is + `test_refine_golden.py`'s CLI quarter circle, a same-CRS domain. +- **Axis order, able to fail:** the run that meshes a 4326 domain unpatched is + refused, naming 4326 and 25833, when `Transformer.from_crs` is patched to + force `always_xy=False`. A GeoJSON written latitude first is refused the same + way. + +**Existing tests changed.** `test_refine_golden.py` (through it, three tests +of `test_cli_constraint_feet.py`) calls `read_domain(path)`; red until green, +for the signature only. In `test_cli_mesh_domain.py` the two tests that met +16's must-match rule are renamed for what now refuses them, the extent check +after the transform, with their assertions unchanged. + ## Test data **Norway (15a, 15b).** diff --git a/docs/increments/16-domain-polygon.md b/docs/increments/16-domain-polygon.md index 601731ef..63007f1b 100644 --- a/docs/increments/16-domain-polygon.md +++ b/docs/increments/16-domain-polygon.md @@ -164,7 +164,9 @@ its corners are NoData vertices. See R5. area-registered DEM this is half a cell inside the image edge. See U4 for clipping instead. The quarter circle's straight edges lie on the border row and column exactly, so they pass. -- **CRS, required, never transformed.** GeoJSON under RFC 7946 is WGS 84 unless +- **CRS, required, never transformed** *(superseded by increment 15b: the + domain is reprojected into the DEM's CRS; see "Ruled by the user", U1).* + GeoJSON under RFC 7946 is WGS 84 unless it carries the 2008 spec's `crs` member, so: a GeoJSON file's CRS is its `crs` member (`urn:ogc:def:crs:EPSG::25833` or `EPSG:25833`), and a file without one is EPSG:4326 by the standard and refused as a mismatch. WKT has no @@ -173,7 +175,8 @@ its corners are NoData vertices. See R5. - **`GeoPolygon` and `--bbox`.** Increment 15 says `GeoPolygon` lands with the clip, as its first caller. This is the clip, but it needs no transform, no `intersects` and no `buffer` yet. So this increment lands a smaller - `DomainPolygon` (shapely polygon plus EPSG code, frozen) in a new + `DomainPolygon` (shapely polygon plus EPSG code, frozen; since 15b, plus a + `crs` string instead of the code) in a new `tin_engine/domain.py`. When 15 lands, its `plan_mosaic` takes `domain.polygon.bounds` in place of `--bbox`; `--domain` and `--bbox` are then mutually exclusive. `GeoPolygon` with `transform` grows from `DomainPolygon` @@ -341,7 +344,11 @@ each (T-deg). - **2026-09-26, the user on U1: "I don't think the domain CRS should have to match the DEM CRS in the future. This must be written down. The questions will be in which coordinate system we shall do the math."** So U1 (a)'s - must-match rule is this increment's scope, not a design principle. A later + must-match rule is this increment's scope, not a design principle. + **Superseded by increment 15b** (`docs/increments/15-dem-mosaic.md`, R9): + the domain is read in its own CRS and reprojected into the DEM's (vertices + only, one `from_crs` site in `crs.py`); for a projected DEM the DEM's CRS is + the computation CRS. A later increment lets the domain, feature geometry and DEM each come in their own CRS, and must first rule on the **computation CRS**: the one coordinate system the noder, CDT, refinement and predicates work in, which today is the diff --git a/project_structure.md b/project_structure.md index a9008825..bfd29b9c 100644 --- a/project_structure.md +++ b/project_structure.md @@ -85,9 +85,12 @@ src_python/tin_engine/ # public Python API (distribution name: rasputin) # pure numpy, never imports _core mosaic.py # plan_mosaic / assemble: select, group by lattice, # check overlaps and coverage, stitch; no files (15a) - dem_input.py # --dem/--bbox -> DemInput(tile, plan, label) (15a) - domain.py # --domain: reads one polygon (GeoJSON or WKT), checks - # CRS and extent; shapely + pyproj, never imports _core + dem_input.py # --dem/--bbox or a domain -> DemInput(tile, plan, + # label, domain in the DEM's CRS) (15a, 15b) + domain.py # --domain: reads one polygon (GeoJSON or WKT) in its + # own CRS, to_crs, check_extent; never imports _core + crs.py # parse_crs, reprojector: the one Transformer.from_crs + # site, always_xy (15b); pyproj and numpy elevation.py # drops mesh vertices the DEM has no data for; # pure numpy, never imports _core stats.py # --stats: PhaseClock, quality, Report, render to @@ -243,9 +246,11 @@ optional NoData sentinel — nothing else. separate from the one above: that one is about where a string may live, this one is about what the numbers mean. Python must deliver every input in one projected, metre CRS, and rejects a geographic CRS outright. Today it does - this by refusal alone: `io/geotiff.py` refuses any file that is not already - in one, and nothing reprojects yet. Reprojection arrives no earlier than - the mosaic increment (`docs/increments/11-raster-ingestion.md` §9, §10). + this partly by refusal: `io/geotiff.py` refuses a DEM that is not already in + one, and since increment 15b the domain polygon is reprojected into the DEM's + CRS in Python (`crs.py`, the one `from_crs` site). A geographic DEM's + computation frame is 15c (`docs/increments/15-dem-mosaic.md`; + `docs/increments/11-raster-ingestion.md` §9, §10). Measured, pyproj 3.8.0 / PROJ 9.8.1: a 0.0002777° cell at 60°N is 15.5 m east-west and 31.0 m north-south, a 2:1 anisotropy invisible to `sample.hpp`, whose bilinear weights would then be computed in degrees and applied to metres. The legacy diff --git a/src_python/tin_engine/cli.py b/src_python/tin_engine/cli.py index 34fb2b31..69b5030e 100644 --- a/src_python/tin_engine/cli.py +++ b/src_python/tin_engine/cli.py @@ -68,6 +68,7 @@ sample, triangulate, ) +from tin_engine.crs import crs_label, parse_crs, transform_description from tin_engine.dem_input import DemInput, DemRequest, open_dem from tin_engine.domain import DomainError, DomainPolygon, read_domain from tin_engine.elevation import Trimmed, trim @@ -596,7 +597,10 @@ def mesh( ] = None, domain_crs: Annotated[ str | None, - typer.Option("--domain-crs", help="The --domain file's CRS, EPSG:n; required for .wkt."), + typer.Option( + "--domain-crs", + help="The --domain file's CRS, anything pyproj reads; required for .wkt.", + ), ] = None, start_min_angle: Annotated[ float | None, @@ -674,6 +678,8 @@ def mesh( raise typer.BadParameter("applies only with --dem", param_hint="--bbox") if domain is None and domain_crs is not None: raise typer.BadParameter("applies only with --domain", param_hint="--domain-crs") + if domain is not None and bbox is not None: + raise typer.BadParameter("--bbox and --domain exclude each other", param_hint="--bbox") if name is not None and name not in GALLERY: raise typer.BadParameter(f"unknown fixture {name}; the gallery is: {', '.join(GALLERY)}") if out.suffix not in MESH_SUFFIXES: @@ -720,7 +726,14 @@ def mesh( raise typer.BadParameter( "--no-constraint-feet needs --tolerance", param_hint="--no-constraint-feet" ) - opened = _open_dem(dem, bbox, clock) + given = None + if domain is not None: + with clock.phase("domain read"): + try: + given = read_domain(domain, domain_crs) + except DomainError as exc: + raise typer.BadParameter(str(exc), param_hint="--domain") from exc + opened = _open_dem(dem, bbox, clock, given) label = opened.label dem_run = _dem_mesh( opened.tile, @@ -730,8 +743,8 @@ def mesh( snap_spacing, tolerance, clock, - domain, - domain_crs, + opened.domain, + domain.name if domain is not None else "", DEFAULT_START_MIN_ANGLE if start_min_angle is None else start_min_angle, not no_constraint_feet, ) @@ -750,9 +763,12 @@ def mesh( seams = opened.seams listed = "; ".join(s.entry() for s in seams) or "none" fields.append(("dem_seams", _ascii(listed))) # Ola's Q1 revised - if described: + if described and given is not None: fields.append(("domain", described)) comments.append(f"domain {described}") + same = parse_crs(given.crs) == parse_crs(f"EPSG:{epsg}") + how = "none" if same else transform_description(given.crs, f"EPSG:{epsg}") + fields += [("domain_crs", crs_label(given.crs)), ("domain_transform", how)] else: assert name is not None if stride is not None: @@ -921,10 +937,14 @@ def _fixture_mesh(name: str, delaunay: bool, spacing: float, clock: PhaseClock) def _open_dem( - dem: list[Path], bbox: tuple[float, float, float, float] | None, clock: PhaseClock + dem: list[Path], + bbox: tuple[float, float, float, float] | None, + clock: PhaseClock, + domain: DomainPolygon | None = None, ) -> DemInput: - """``--dem`` and ``--bbox`` to one tile (increment 15a, R11); every refusal, - the reader's or the mosaic's, is a usage error in its own words.""" + """``--dem`` and ``--bbox`` or the read ``--domain`` to one tile (increment + 15a and 15b, R11); every refusal, the reader's, the mosaic's or the + domain's extent, is a usage error in its own words.""" try: bounds = ( None @@ -935,12 +955,14 @@ def _open_dem( raise typer.BadParameter(_words(exc), param_hint="--bbox") from exc try: with clock.phase("decode"): - return open_dem(DemRequest(sources=tuple(dem), bounds=bounds)) + return open_dem(DemRequest(sources=tuple(dem), bounds=bounds, domain=domain)) except OSError as exc: where = exc.filename or ", ".join(map(str, dem)) raise typer.BadParameter( f"cannot read {where}: {exc.strerror or exc}", param_hint="--dem" ) from exc + except DomainError as exc: + raise typer.BadParameter(str(exc), param_hint="--domain") from exc except ValueError as exc: raise typer.BadParameter(_words(exc), param_hint="--dem") from exc @@ -982,8 +1004,8 @@ def _dem_mesh( spacing: float, tolerance: float | None, clock: PhaseClock, - domain: Path | None = None, - domain_crs: str | None = None, + domain: DomainPolygon | None = None, + domain_name: str = "", min_angle: float = 0.0, feet: bool = False, ) -> _DemMesh: @@ -991,8 +1013,9 @@ def _dem_mesh( Without ``tolerance`` this is increment 12's R6: z sampled bilinearly at the stride grid. With it, increment 14's R9: the stride grid is the start - mesh, refined against the DEM's nodes; with ``domain``, increment 16's R3, - the polygon's rings are. ``min_angle`` > 0 improves the start's angles + mesh, refined against the DEM's nodes; with ``domain`` (already in the + DEM's CRS, 15b; ``domain_name`` is its file's), increment 16's R3, the + polygon's rings are. ``min_angle`` > 0 improves the start's angles first (increment 20); ``feet`` inserts constraint feet (increment 20b). Returns the mesh, the ``elevation`` sentence for the file, the EPSG code, the ``domain`` field (empty without one), and the ``--stats`` inputs; @@ -1004,13 +1027,8 @@ def _dem_mesh( described = "" domain_vertices = domain_holes = None if domain is not None: - with clock.phase("domain read"): - try: - polygon = read_domain(domain, meta, domain_crs) - except DomainError as exc: - raise typer.BadParameter(str(exc), param_hint="--domain") from exc - xy, chains, described = _domain_chains(polygon, domain.name) - domain_vertices, domain_holes = len(xy), len(polygon.polygon.interiors) + xy, chains, described = _domain_chains(domain, domain_name) + domain_vertices, domain_holes = len(xy), len(domain.polygon.interiors) start = "start domain boundary, boundary z bilinear" run = _engine(xy, chains, delaunay, spacing, clock) else: diff --git a/src_python/tin_engine/crs.py b/src_python/tin_engine/crs.py new file mode 100644 index 00000000..536a717c --- /dev/null +++ b/src_python/tin_engine/crs.py @@ -0,0 +1,67 @@ +"""CRS text and the one reprojection site (increment 15b). + +`docs/increments/15-dem-mosaic.md` R2, R9 and I8. `_transformer` holds the only +`Transformer.from_crs` in `src_python/`, always with `always_xy=True`: x is +easting or longitude and y northing or latitude, whatever the CRS's own axis +order, which is the convention a GeoTIFF's model space and GeoJSON both use. + +Every refusal is a `ValueError`. pyproj's `CRSError` is a `RuntimeError`, and +would escape every `except ValueError` that turns a refusal into a usage error. + +Pure: pyproj and numpy. No paths, no `_core`; CRS never crosses into the core. +""" + +from __future__ import annotations + +from collections.abc import Callable +from typing import Any + +import numpy as np +import numpy.typing as npt +from pyproj import CRS, Transformer +from pyproj.exceptions import CRSError + +Xy = npt.NDArray[np.float64] + + +def parse_crs(text: str | CRS) -> CRS: + """Anything `CRS.from_user_input` accepts, or a `ValueError` naming it.""" + if isinstance(text, CRS): + return text + try: + return CRS.from_user_input(text) + except CRSError as exc: + raise ValueError(f"cannot read the CRS {text!r}: {exc}") from exc + + +def reprojector(src: str | CRS, dst: str | CRS) -> Callable[[Any], Xy]: + """A map from an `(N, 2)` array-like of `(x, y)` in `src` to an `(N, 2)` + float64 array in `dst`: pyproj's `always_xy` transform, bit for bit. + A point with no image comes back as `inf`; the caller decides.""" + transformer = _transformer(src, dst) + + def apply(xy: Any) -> Xy: + points = np.asarray(xy, dtype=np.float64).reshape(-1, 2) + x, y = transformer.transform(points[:, 0], points[:, 1]) + return np.column_stack([np.asarray(x, np.float64), np.asarray(y, np.float64)]) + + return apply + + +def transform_description(src: str | CRS, dst: str | CRS) -> str: + """What PROJ picked for `src` to `dst` (a datum shift included), for the record.""" + return str(_transformer(src, dst).description) + + +def crs_label(crs: str | CRS) -> str: + """`EPSG:n` when pyproj finds an exact EPSG code (so OGC's CRS84 is not + `EPSG:4326`), otherwise ASCII text pyproj parses back to the same CRS.""" + parsed = parse_crs(crs) + code = parsed.to_epsg(min_confidence=100) + if code is not None: + return f"EPSG:{code}" + return parsed.to_string().encode("ascii", "backslashreplace").decode("ascii") + + +def _transformer(src: str | CRS, dst: str | CRS) -> Transformer: + return Transformer.from_crs(parse_crs(src), parse_crs(dst), always_xy=True) diff --git a/src_python/tin_engine/dem_input.py b/src_python/tin_engine/dem_input.py index 43a8cd37..01895384 100644 --- a/src_python/tin_engine/dem_input.py +++ b/src_python/tin_engine/dem_input.py @@ -6,29 +6,38 @@ `open_dem`; a GUI backend or an API worker builds the same request without Typer. Besides `io/repository.py`, this is the one module below `cli.py` with paths, and it only hands them to the repository. + +With a domain (increment 15b, R4 point 5, R6 and R9), the domain is moved into +the DEM's CRS first; its bounds take the box's place, the polygon grown by one +cell is the region whose nodes must be covered, and its extent is checked +against the plan, all before any tile is loaded (I6). """ from __future__ import annotations +import math from dataclasses import dataclass from pathlib import Path -from typing import Self +from typing import Any, Self from pydantic import BaseModel, ConfigDict, model_validator -from tin_engine.io.models import DemTile +from tin_engine.domain import DomainError, DomainPolygon, check_extent +from tin_engine.io.models import DemTile, RasterMeta from tin_engine.io.repository import TiffDemRepository -from tin_engine.mosaic import Bounds, MosaicPlan, Seam, assemble, plan_mosaic +from tin_engine.mosaic import Bounds, MosaicError, MosaicPlan, Seam, assemble, plan_mosaic class DemRequest(BaseModel): - """Exactly one directory of tiles, or one or more tile files (R11).""" + """Exactly one directory of tiles, or one or more tile files (R11), and at + most one of a box in the DEM's CRS and a domain in its own (R6).""" model_config = ConfigDict(frozen=True) sources: tuple[Path, ...] bounds: Bounds | None = None nodata: float | None = None + domain: DomainPolygon | None = None @model_validator(mode="after") def _one_form(self) -> Self: @@ -40,18 +49,22 @@ def _one_form(self) -> Self: f"give exactly one directory, or one or more files; got {len(self.sources)} " f"sources of which {directories} are directories" ) + if self.bounds is not None and self.domain is not None: + raise ValueError("give a box or a domain, not both: the domain's bounds are the box") return self @dataclass(frozen=True, slots=True) class DemInput: """The tile to mesh, the plan it was assembled by, a name for the run (the - directory's name, or the stem of the first file as given), and the - mosaic's disagreeing seams (Ola's Q1 revised).""" + directory's name, or the stem of the first file as given), the domain in + the DEM's CRS (None without one), and the mosaic's disagreeing seams (Ola's + Q1 revised).""" tile: DemTile plan: MosaicPlan label: str + domain: DomainPolygon | None = None seams: tuple[Seam, ...] = () @@ -65,6 +78,61 @@ def open_dem(request: DemRequest) -> DemInput: else: repository = TiffDemRepository(request.sources, nodata=request.nodata) label = first.stem - plan = plan_mosaic(repository.footprints(), request.bounds, None) - mosaic = assemble(plan, repository.load) - return DemInput(tile=mosaic.tile, plan=plan, label=label, seams=mosaic.seams) + footprints = repository.footprints() + if request.domain is None or not footprints: + plan, domain, grown = plan_mosaic(footprints, request.bounds, None), None, None + else: + plan, domain, grown = _domain_plan(footprints, request.domain) + # With a domain the seam report counts only the needed region's nodes + # (Ola, 2026-09-28); which value a node takes does not depend on it. + mosaic = assemble(plan, repository.load, grown) + return DemInput(tile=mosaic.tile, plan=plan, label=label, domain=domain, seams=mosaic.seams) + + +def _domain_plan(footprints: Any, given: DomainPolygon) -> tuple[MosaicPlan, DomainPolygon, Any]: + """The domain in the DEM's CRS, the plan of its bounds with the polygon + grown by one cell needed (R4 point 5), and that grown polygon, the needed + region. Grown by the cell diagonal with mitred corners, a superset of the + four nodes bilinear z reads at any point of the polygon. A refusal names + both CRSs and keeps its type.""" + epsgs = sorted({f.meta.epsg for f in footprints}) + if len(epsgs) > 1: + raise MosaicError(f"the tiles are in {len(epsgs)} CRSs, EPSG:{epsgs}; a domain needs one") + try: + domain = given.to_crs(f"EPSG:{epsgs[0]}") + box = Bounds( + **dict(zip(("x_min", "y_min", "x_max", "y_max"), domain.polygon.bounds, strict=True)) + ) + # The cell is the chosen lattice's, not a listed tile's (15b review, B1): + # plan on the polygon itself, grow by that plan's cell diagonal, and + # re-plan while the lattice chosen has a larger one. The reach only + # grows, so this ends, and the region needed is never short of the + # chosen lattice's own cell. + plan, reach, grown = plan_mosaic(footprints, box, domain.polygon), 0.0, domain.polygon + while (diagonal := math.hypot(plan.meta.delta_x, plan.meta.delta_y)) > reach: + reach, grown = diagonal, domain.polygon.buffer(diagonal, join_style="mitre") + plan = plan_mosaic(footprints, box, grown) + if past := _past(box, plan.meta): # a vertex within 1e-6 cell past a node line + plan = plan_mosaic(footprints, past, grown) + check_extent(domain, plan.meta) + except (DomainError, MosaicError) as exc: + raise type(exc)(f"the domain, in {given.crs}, in the DEM's EPSG:{epsgs[0]}: {exc}") from exc + return plan, domain, grown + + +def _past(box: Bounds, m: RasterMeta) -> Bounds | None: + """`box` with each edge past `m`'s node rectangle moved out by half a cell, + or None when none is. 15a's window snaps an edge within `ALIGN_TOLERANCE` + cell of a node line onto it; a domain vertex just past that line needs the + next one, and moving the edge half a cell out takes exactly that line.""" + x_max, y_min = m.x_min + (m.cols - 1) * m.delta_x, m.y_max - (m.rows - 1) * m.delta_y + out = (box.x_min < m.x_min, box.y_min < y_min, box.x_max > x_max, box.y_max > m.y_max) + if not any(out): + return None + hx, hy = m.delta_x / 2, m.delta_y / 2 + return Bounds( + x_min=box.x_min - hx * out[0], + y_min=box.y_min - hy * out[1], + x_max=box.x_max + hx * out[2], + y_max=box.y_max + hy * out[3], + ) diff --git a/src_python/tin_engine/domain.py b/src_python/tin_engine/domain.py index 7831fccc..c5205613 100644 --- a/src_python/tin_engine/domain.py +++ b/src_python/tin_engine/domain.py @@ -5,15 +5,17 @@ given: no densifying, simplifying or snapping. The polygon is oriented outer counter-clockwise and holes clockwise (increment 3's winding contract). -CRS is required and never transformed. A GeoJSON file's is its ``crs`` member, -and a file without one is EPSG:4326 by RFC 7946; WKT has none, so it comes from -the caller. The must-match rule is this increment's scope only (U1, and the -user's caveat on it), so it lives in one replaceable function, :func:`check_crs`. +CRS is required. A GeoJSON file's is its ``crs`` member, and a file without one +is EPSG:4326 by RFC 7946; WKT has none, so it comes from the caller. The domain +keeps its own CRS, and :meth:`DomainPolygon.to_crs` moves it into the DEM's +(increment 15b, ``15-dem-mosaic.md`` R9), replacing 16's must-match rule. Every vertex must lie in the DEM's node rectangle, border included (U4 a): -outside it there is no bilinear z (R0). +outside it there is no bilinear z (R0). That is :func:`check_extent`, run after +the transform, against the mosaic. -Pure: json, shapely, pyproj and ``RasterMeta``. No ``_core``, no typer. +Pure: json, numpy, shapely, pyproj, pydantic, ``tin_engine.crs`` and +``RasterMeta``. No ``_core``, no typer. """ from __future__ import annotations @@ -22,20 +24,21 @@ from pathlib import Path from typing import Any +import numpy as np import shapely import shapely.wkt from pydantic import BaseModel, ConfigDict from pyproj import CRS -from pyproj.exceptions import CRSError from shapely.geometry import Polygon, shape from shapely.geometry.polygon import orient from shapely.validation import explain_validity +from tin_engine.crs import parse_crs, reprojector from tin_engine.io.models import RasterMeta GEOJSON_SUFFIXES = (".geojson", ".json") WKT_SUFFIXES = (".wkt",) -GEOJSON_DEFAULT_EPSG = 4326 # RFC 7946: no crs member means WGS 84 +GEOJSON_DEFAULT_CRS = "EPSG:4326" # RFC 7946: no crs member means WGS 84 class DomainError(ValueError): @@ -43,16 +46,38 @@ class DomainError(ValueError): class DomainPolygon(BaseModel): - """The domain as read: an oriented shapely polygon and its EPSG code.""" + """The domain: an oriented shapely polygon, and the text of its CRS.""" model_config = ConfigDict(frozen=True, arbitrary_types_allowed=True) polygon: Polygon - epsg: int + crs: str + + def to_crs(self, dst: str | CRS) -> DomainPolygon: + """The domain in ``dst``: each vertex through pyproj's ``always_xy`` + transform, once, and re-oriented; rings stay straight (R9). An equal + CRS, by ``pyproj.CRS`` equality and not by EPSG code, returns ``self``. + """ + source, target = parse_crs(self.crs), parse_crs(dst) + if source == target: + return self + move = reprojector(source, target) + rings = [np.asarray(r.coords) for r in (self.polygon.exterior, *self.polygon.interiors)] + moved = [move(r) for r in rings] + for ring, image in zip(rings, moved, strict=True): + bad = ~np.isfinite(image).all(axis=1) + if bad.any(): + x, y = ring[bad][0] + raise DomainError( + f"the domain vertex ({x}, {y}) in {self.crs} has no image in " + f"{target.to_string()}" + ) + polygon = orient(Polygon(moved[0], moved[1:]), sign=1.0) + return DomainPolygon(polygon=polygon, crs=target.to_string()) -def read_domain(path: Path, meta: RasterMeta, crs: str | None = None) -> DomainPolygon: - """Read, check and orient the domain in ``path`` against the DEM ``meta``. +def read_domain(path: Path, crs: str | None = None) -> DomainPolygon: + """Read, check and orient the domain in ``path``, in its own CRS. ``crs`` is the ``--domain-crs`` text: required for WKT, and for GeoJSON it must agree with the file's own. Raises :class:`DomainError` on any refusal. @@ -71,13 +96,12 @@ def read_domain(path: Path, meta: RasterMeta, crs: str | None = None) -> DomainP if suffix in WKT_SUFFIXES: if crs is None: raise DomainError(f"{path.name}: WKT carries no CRS; give it with --domain-crs") - geometry, epsg = shapely.wkt.loads(text), _epsg(crs) + geometry, own = shapely.wkt.loads(text), crs + _parsed(own) else: - geometry, epsg = _geojson(json.loads(text)) - if crs is not None and _epsg(crs) != epsg: - raise DomainError( - f"{path.name} is in EPSG:{epsg} but --domain-crs says EPSG:{_epsg(crs)}" - ) + geometry, own = _geojson(json.loads(text)) + if crs is not None and _parsed(crs) != _parsed(own): + raise DomainError(f"{path.name} is in {own} but --domain-crs says {crs}") except DomainError: raise except (ValueError, TypeError, KeyError, AttributeError, shapely.errors.GEOSException) as exc: @@ -89,35 +113,22 @@ def read_domain(path: Path, meta: RasterMeta, crs: str | None = None) -> DomainP raise DomainError(f"{path.name}: the domain polygon is empty") if not geometry.is_valid: raise DomainError(f"{path.name}: invalid polygon, {explain_validity(geometry)}") - check_crs(epsg, meta) - _check_extent(geometry, meta) - return DomainPolygon(polygon=orient(geometry, sign=1.0), epsg=epsg) - - -def check_crs(epsg: int, meta: RasterMeta) -> None: - """U1 (a): the domain's CRS must be the DEM's, since nothing reprojects yet.""" - if epsg != meta.epsg: - raise DomainError( - f"the domain is in EPSG:{epsg} and the DEM in EPSG:{meta.epsg}; they must match " - f"(a GeoJSON file without a crs member is EPSG:{GEOJSON_DEFAULT_EPSG})" - ) + return DomainPolygon(polygon=orient(geometry, sign=1.0), crs=own) -def _epsg(text: str) -> int: +def _parsed(text: str) -> CRS: try: - code = CRS.from_user_input(text).to_epsg() - except CRSError as exc: - raise DomainError(f"cannot read the CRS {text!r}: {exc}") from exc - if code is None: - raise DomainError(f"the CRS {text!r} has no EPSG code") - return int(code) + return parse_crs(text) + except ValueError as exc: + raise DomainError(str(exc)) from exc -def _geojson(doc: dict[str, Any]) -> tuple[Any, int]: +def _geojson(doc: dict[str, Any]) -> tuple[Any, str]: """The geometry of a bare geometry, a Feature or a one-feature collection, - and the EPSG code of the document's ``crs`` member.""" + and the CRS text of the document's ``crs`` member, checked.""" member = doc.get("crs") - epsg = _epsg(member["properties"]["name"]) if member is not None else GEOJSON_DEFAULT_EPSG + own = str(member["properties"]["name"]) if member is not None else GEOJSON_DEFAULT_CRS + _parsed(own) if doc.get("type") == "FeatureCollection": features = doc["features"] if len(features) != 1: @@ -125,11 +136,13 @@ def _geojson(doc: dict[str, Any]) -> tuple[Any, int]: doc = features[0] if doc.get("type") == "Feature": doc = doc["geometry"] - return shape(doc), epsg + return shape(doc), own -def _check_extent(polygon: Polygon, meta: RasterMeta) -> None: - """U4 (a): every vertex in the node rectangle, as the core's ``cell_of`` has it.""" +def check_extent(domain: DomainPolygon, meta: RasterMeta) -> None: + """U4 (a): every vertex in ``meta``'s node rectangle, as the core's + ``cell_of`` has it, with ``domain`` already in the DEM's CRS.""" + polygon = domain.polygon x_max = meta.x_min + (meta.cols - 1) * meta.delta_x y_min = meta.y_max - (meta.rows - 1) * meta.delta_y for ring in (polygon.exterior, *polygon.interiors): diff --git a/src_python/tin_engine/mosaic.py b/src_python/tin_engine/mosaic.py index ff8d46b2..ee66fd01 100644 --- a/src_python/tin_engine/mosaic.py +++ b/src_python/tin_engine/mosaic.py @@ -226,7 +226,7 @@ def plan_mosaic( ) -def assemble(plan: MosaicPlan, load: Callable[[str], DemTile]) -> Mosaic: +def assemble(plan: MosaicPlan, load: Callable[[str], DemTile], needed: Any = None) -> Mosaic: """Load the plan's tiles one at a time, in plan order, into one canvas (R5). One tile whose grid is the mosaic's is returned as loaded: no canvas, no @@ -240,6 +240,11 @@ def assemble(plan: MosaicPlan, load: Callable[[str], DemTile]) -> Mosaic: are in, every overlap is decided from the strips (Ola's Q1 revised, `_decide`) and each pair that disagrees is reported (`_seam`). + `needed`, a shapely geometry in the DEM's CRS as for `plan_mosaic`, limits + the report to the overlap nodes it covers (closed; Ola, 2026-09-28: with + `--domain`, the needed region). Only the report: every node is decided as + without it, and only overlap nodes are tested, never the canvas. + Memory: the peak is the canvas, the strips, and one load's own peak, which is not one tile: decoding a DTM10 tile peaks at 2.0 to 3.0 tiles, depending on the tile (the decoder's buffers and `DemTile`'s read-only @@ -268,13 +273,17 @@ def assemble(plan: MosaicPlan, load: Callable[[str], DemTile]) -> Mosaic: strips[placement.name, other.name] = incoming[_within(box, c)].copy() del array, incoming # before the next load, or two tiles outlive this one seams = [] + if needed is not None: + shapely.prepare(needed) for i, a in enumerate(ordered): for b in ordered[i + 1 :]: box = _meet(a.canvas, b.canvas) if box is None: continue _decide(canvas, box, ordered, strips, (a.name, b.name), plan.meta.nodata) - seam = _seam(a.name, b.name, strips[a.name, b.name], strips[b.name, a.name], plan) + kept = _covered(box, plan.meta, needed) + first, second = strips[a.name, b.name][kept], strips[b.name, a.name][kept] + seam = _seam(a.name, b.name, first, second, plan) if seam is not None: seams.append(seam) return Mosaic(tile=DemTile._adopt(plan.meta, canvas), plan=plan, seams=tuple(seams)) @@ -519,6 +528,18 @@ def _decide( canvas[box.row0 : box.row0 + box.rows, box.col0 : box.col0 + box.cols] = value +def _covered(box: IndexWindow, meta: RasterMeta, needed: Any) -> Any: + """The nodes of `box` that `needed` covers, as a boolean mask of its shape, + or `...` (every node) without it. A node is `x_min + col * dx`, as in + `_uncovered`.""" + if needed is None: + return ... + rows = (box.row0 + np.arange(box.rows))[:, np.newaxis] + cols = (box.col0 + np.arange(box.cols))[np.newaxis, :] + xs, ys = np.broadcast_arrays(meta.x_min + cols * meta.delta_x, meta.y_max - rows * meta.delta_y) + return shapely.intersects_xy(needed, xs, ys) + + def _seam( first: str, second: str, a: npt.NDArray[Any], b: npt.NDArray[Any], plan: MosaicPlan ) -> Seam | None: diff --git a/tests/python/mosaic_fixtures.py b/tests/python/mosaic_fixtures.py index ee6c9539..310efd1b 100644 --- a/tests/python/mosaic_fixtures.py +++ b/tests/python/mosaic_fixtures.py @@ -234,20 +234,25 @@ def winners(tiles: Mapping[str, DemTile], grid: RasterMeta) -> list[list[str]]: def seams_of( - tiles: Mapping[str, DemTile], grid: RasterMeta + tiles: Mapping[str, DemTile], grid: RasterMeta, needed: np.ndarray | None = None ) -> list[tuple[str, str, int, float, float]]: """The seam report the rule asks for, pair by pair, sorted by name: `(first, second, nodes, largest, median)` over the mosaic's nodes where both tiles hold a valid value and `|a - b| >= SEAM_THRESHOLD`; `|difference|` in float64, the median of an even count the mean of the - middle two. Pairs with no such node are left out.""" + middle two. Pairs with no such node are left out. + + `needed`, a boolean array of `grid`'s shape, keeps only the nodes it marks + (Ola, 2026-09-28: with `--domain`, the nodes of the needed region).""" placed = on_canvas(tiles, grid) + kept = np.ones((grid.rows, grid.cols), dtype=bool) if needed is None else needed + assert kept.shape == (grid.rows, grid.cols), (kept.shape, grid.rows, grid.cols) out = [] names = sorted(placed) for i, a in enumerate(names): for b in names[i + 1 :]: va, vb = placed[a][0], placed[b][0] - both = ~np.isnan(va) & ~np.isnan(vb) + both = kept & ~np.isnan(va) & ~np.isnan(vb) differ = both & (np.abs(np.where(both, va - vb, 0.0)) >= SEAM_THRESHOLD) if differ.any(): gaps = np.abs(va[differ] - vb[differ]) diff --git a/tests/python/test_cli_mesh_domain.py b/tests/python/test_cli_mesh_domain.py index 17fff5a7..602ea222 100644 --- a/tests/python/test_cli_mesh_domain.py +++ b/tests/python/test_cli_mesh_domain.py @@ -2,9 +2,13 @@ `docs/increments/16-domain-polygon.md`, R1 to R4, with the user's rulings R0 (point heights, bilinear between nodes), no snapping of input vertices, U1 (a) -(the domain's CRS must match the DEM's), U3 (a), U4 (a) (refuse a polygon -outside the node rectangle), U5 (a) (no boundary bit) and U6 (a) (the noder's -1 mm snap is the engine's input precision). +(the domain's CRS is its `crs` member, or `--domain-crs` for WKT), U3 (a), +U4 (a) (refuse a polygon outside the node rectangle), U5 (a) (no boundary bit) +and U6 (a) (the noder's 1 mm snap is the engine's input precision). U1 (a)'s +must-match rule is replaced by increment 15b's transform +(`15-dem-mosaic.md` R9; `test_cli_mesh_domain_crs.py`): the two tests below +that met it now meet the extent check after the transform, which names both +CRSs. Wording pinned from the design: the ``domain`` field is ``, 1 ring holes, vertices`` (``hole`` or ``holes`` @@ -353,11 +357,13 @@ def test_empty(self, tmp_path: Path, bumpy: Path) -> None: says=("--domain",), ) - def test_geojson_without_crs(self, tmp_path: Path, bumpy: Path) -> None: + def test_utm_numbers_without_a_crs_member(self, tmp_path: Path, bumpy: Path) -> None: + """WGS 84 by RFC 7946, so "longitude" 500 012: no image in UTM 33 (15b).""" path = geojson(tmp_path / "wgs.geojson", SQUARE, crs=None) self.refused(tmp_path, bumpy, path, "--tolerance", "1", says=("4326", "25833")) - def test_a_mismatched_epsg(self, tmp_path: Path, bumpy: Path) -> None: + def test_utm33_numbers_labelled_utm32_land_outside(self, tmp_path: Path, bumpy: Path) -> None: + """Transformed from zone 32 (15b), they land about 340 km west of the DEM.""" path = geojson(tmp_path / "utm32.geojson", SQUARE, crs="EPSG:25832") self.refused(tmp_path, bumpy, path, "--tolerance", "1", says=("25832", "25833")) diff --git a/tests/python/test_cli_mesh_domain_crs.py b/tests/python/test_cli_mesh_domain_crs.py new file mode 100644 index 00000000..08319cee --- /dev/null +++ b/tests/python/test_cli_mesh_domain_crs.py @@ -0,0 +1,495 @@ +"""`rasputin mesh --dem ... --domain PATH` with the domain in its own CRS (increment 15b). + +`docs/increments/15-dem-mosaic.md` R9, R11 and "Tests for @tester" (15b), and +the 15b Acceptance: a catchment polygon in EPSG:4326 over a 25833 mosaic +meshes, with `domain_crs` and `domain_transform` recorded. Pinned by this suite +(see "Pinned by the red suite (15b)"): + +- `domain_crs` is `EPSG:n` when pyproj finds an exact EPSG code for the + domain's CRS, and otherwise text pyproj parses back to that CRS. +- `domain_transform` is pyproj's `description` of the `always_xy` transformer + from the domain's CRS to the DEM's; for a domain already in the DEM's CRS it + says `none` (any case), and no transformer is made. +- A domain already in the DEM's CRS meshes bit for bit as increment 16 meshed + it, and no transformer is made. Relational, on the same machine: the run + equals, in every point, cell and array, the same run with `cli.open_dem` + replaced by 16's data flow (the whole file, and the domain object exactly as + `read_domain` returned it, never through `DomainPolygon.to_crs`), and the + domain `_dem_mesh` receives has the read vertices bit for bit. This + replaced two digests recorded on macOS arm64 at `d34d79d` (test amendment + after PR #106's CI): Linux x86 with GCC gives the square a different digest, + so a recorded digest pins a platform, not a behaviour. The platform-stable + absolute anchor for a same-CRS domain through the CLI is increment 18's + `test_refine_golden.py::test_the_cli_with_start_min_angle_0_matches_the_digest` + (the quarter circle, whose GeoJSON names EPSG:25833); it is referenced, not + duplicated. +- `--bbox` with `--domain` is a usage error naming both. +- `dem_seams` is written on the domain path as without a domain: the + disagreeing pair with a domain in EPSG:4326, `none` when the overlaps agree + (test amendment after the 15b review). With a domain it counts only the + nodes inside the needed region (Ola, 2026-09-28): `none` for a disagreement + wholly outside it, however much of the plan's rectangle it fills. +- A transformed domain refused for its extent (outside the DEM, in no tile, or + with no image in the DEM's CRS) is a usage error naming the domain's CRS and + the DEM's EPSG code, and writes nothing. + +The axis-order test is able to fail: it patches `pyproj.Transformer.from_crs` +to drop `always_xy`, and the same run that meshes unpatched must then be +refused by the extent check (R9, "Degeneracy policy": axis order). + +HOW THIS FILE GOES RED: it imports nothing new. Before 15b a domain in +another CRS is refused as a CRS mismatch, `domain_crs` is not written, and +`--bbox` with `--domain` is accepted, so those tests fail on their assertions. +Six are guards that pass before 15b and must stay green after it: the two +same-CRS runs (then recorded as digests, now relational), the two unknown +`--domain-crs` refusals, and the far-away and latitude-first refusals (16's +mismatch message already names both CRSs). +`test_axis_order_is_able_to_fail` is red before 15b through its control run. +""" + +from __future__ import annotations + +import dataclasses +import hashlib +import json +from pathlib import Path +from typing import Any, ClassVar + +import numpy as np +import pytest +import shapely +from pyproj import CRS, Transformer +from shapely.geometry import Polygon + +import test_cli_mesh_mosaic +import test_dem_input_domain +import tin_engine.cli as cli +from geotiff_fixtures import KARTVERKET, micro_tiff, needs_codecs +from mosaic_fixtures import X0, Y0, blocks, quadrants, whole +from test_cli_mesh_dem import write_tiff +from test_cli_mesh_domain import COLS, HOLE, ROWS, SQUARE, quarter_circle +from test_cli_mesh_domain import geojson as utm33_geojson +from test_cli_mesh_mosaic import USAGE, invoke, terrain, write_tiles +from test_cli_mesh_refine import NUMBER, sentence +from test_cli_mesh_refine import field as sentence_field +from tin_engine.dem_input import DemInput, DemRequest, open_dem +from tin_engine.domain import DomainPolygon, read_domain +from vtkread import VtkFile, read_vtk + +Ring = list[tuple[float, float]] +SNAP = 1e-3 # DEFAULT_SNAP_SPACING: the noder's grid, applied in the DEM's CRS (R9) +LCC = "+proj=lcc +lat_1=60 +lat_2=65 +lat_0=62 +lon_0=15 +ellps=GRS80 +units=m +no_defs" + + +def digest(vtk: VtkFile) -> str: + """SHA-256 over the mesh itself: points, cells and every cell and point + array, each prefixed with its name, dtype and shape. Field data is left out, + since 15b adds `domain_crs` and `domain_transform` to it.""" + h = hashlib.sha256() + + def put(name: str, a: Any) -> None: + arr = np.ascontiguousarray(np.asarray(a)) + h.update(f"{name}{arr.dtype.str}{arr.shape}".encode()) + h.update(arr.tobytes()) + + put("points", vtk.points) + for i, cell in enumerate(vtk.lines): + put(f"line{i}", cell) + for i, cell in enumerate(vtk.polygons): + put(f"polygon{i}", cell) + for group in (vtk.point_scalars, vtk.scalars): + for name in sorted(group): + put(name, group[name].values) + for block in sorted(vtk.cell_fields): + for name in sorted(vtk.cell_fields[block]): + values = vtk.cell_fields[block][name].values + if isinstance(values, np.ndarray): + put(f"{block}.{name}", values) + return h.hexdigest() + + +def to_crs(src: str, dst: str, ring: Ring) -> Ring: + """pyproj's own `always_xy` transform of `ring`: the oracle.""" + t = Transformer.from_crs(src, dst, always_xy=True) + xs, ys = t.transform(np.array([p[0] for p in ring]), np.array([p[1] for p in ring])) + return [(float(x), float(y)) for x, y in zip(xs, ys, strict=True)] + + +def description(src: str, dst: str = "EPSG:25833") -> str: + return str(Transformer.from_crs(src, dst, always_xy=True).description) + + +def write_geojson( + path: Path, outer: Ring, holes: tuple[Ring, ...] = (), crs: str | None = None +) -> Path: + doc: dict[str, Any] = {"type": "Polygon", "coordinates": [[*r, r[0]] for r in (outer, *holes)]} + if crs is not None: + doc["crs"] = {"type": "name", "properties": {"name": crs}} + path.write_text(json.dumps(doc)) + return path + + +def the_field(vtk: VtkFile, name: str) -> str: + assert name in vtk.field_data, sorted(vtk.field_data) + (value,) = vtk.field_data[name].values + return str(value) + + +def run(tmp_path: Path, *args: str, out: str = "x.vtk") -> VtkFile: + target = tmp_path / out + code, output = invoke(*args, "--out", str(target)) + assert code == 0, output + return read_vtk(target.read_bytes()) + + +def refused(tmp_path: Path, *args: str, says: tuple[str, ...]) -> str: + target = tmp_path / "refused.vtk" + code, output = invoke(*args, "--out", str(target)) + assert code == USAGE, output + assert "Traceback" not in output + assert "No such option" not in output, "refused for the wrong reason" + for word in says: + assert word in output, f"{word!r} not in {output!r}" + assert not target.exists() + return output + + +def assert_ring_in_output(vtk: VtkFile, ring: Ring) -> None: + """Every expected vertex has an output point within the noder's snap (U6).""" + points = vtk.points[:, :2] + for x, y in ring: + i = int(np.argmin(np.hypot(points[:, 0] - x, points[:, 1] - y))) + assert abs(points[i, 0] - x) <= SNAP / 2 + 1e-9, (x, y, points[i]) + assert abs(points[i, 1] - y) <= SNAP / 2 + 1e-9, (x, y, points[i]) + + +def assert_inside(vtk: VtkFile, ring: Ring) -> None: + """No output vertex outside the domain, in the DEM's CRS, beyond the snap.""" + grown = Polygon(ring).buffer(SNAP) + inside = shapely.intersects_xy(grown, vtk.points[:, 0], vtk.points[:, 1]) + assert inside.all(), vtk.points[~inside][:5] + + +# ------------------------------------------------------------------ tiles + +# `whole(9, 13)` of `terrain`, in quadrants: nodes x 500 000 .. 500 120 (dx 10), +# y 6 600 000 .. 6 599 960 (dy 5), EPSG:25833, point-registered. +ACROSS: Ring = [ + (X0 + 12.3, Y0 - 36.3), + (X0 + 107.7, Y0 - 35.9), + (X0 + 106.1, Y0 - 3.7), + (X0 + 13.9, Y0 - 4.1), +] + + +@pytest.fixture +def quad_dir(tmp_path: Path) -> Path: + tiles = quadrants(whole(9, 13, array=terrain(9, 13)), row_cut=4, col_cut=6, overlap=1) + write_tiles(tmp_path / "quad", tiles) + return tmp_path / "quad" + + +def mesh_args(dem: Path, domain: Path, *extra: str) -> tuple[str, ...]: + return ("--dem", str(dem), "--domain", str(domain), "--tolerance", "1", *extra) + + +# ------------------------------------------------------------------ tests + + +class TestTheDomainInItsOwnCrs: + """ "Tests for @tester", 15b: 4326 GeoJSON, 25832, WKT with --domain-crs + EPSG:3035, each over a 25833 mosaic; plus a CRS with no EPSG code and + OGC's CRS84.""" + + def test_wgs84_geojson_without_a_crs_member(self, tmp_path: Path, quad_dir: Path) -> None: + lon_lat = to_crs("EPSG:25833", "EPSG:4326", ACROSS) + domain = write_geojson(tmp_path / "catchment.geojson", lon_lat) + vtk = run(tmp_path, *mesh_args(quad_dir, domain)) + expected = to_crs("EPSG:4326", "EPSG:25833", lon_lat) + assert_ring_in_output(vtk, expected) + assert_inside(vtk, expected) + assert the_field(vtk, "domain_crs") == "EPSG:4326" + assert the_field(vtk, "domain_transform") == description("EPSG:4326") + assert the_field(vtk, "crs") == "EPSG:25833" + assert the_field(vtk, "dem_tiles") == "ne.tif; nw.tif; se.tif; sw.tif" + assert sentence_field(sentence(vtk), rf"achieved max error {NUMBER} m") <= 1.0 + + def test_utm32_geojson(self, tmp_path: Path, quad_dir: Path) -> None: + ring = to_crs("EPSG:25833", "EPSG:25832", ACROSS) + domain = write_geojson(tmp_path / "d.geojson", ring, crs="urn:ogc:def:crs:EPSG::25832") + vtk = run(tmp_path, *mesh_args(quad_dir, domain)) + assert_ring_in_output(vtk, to_crs("EPSG:25832", "EPSG:25833", ring)) + assert the_field(vtk, "domain_crs") == "EPSG:25832" + assert the_field(vtk, "domain_transform") == description("EPSG:25832") + + def test_wkt_in_laea_europe_with_domain_crs(self, tmp_path: Path, quad_dir: Path) -> None: + ring = to_crs("EPSG:25833", "EPSG:3035", ACROSS) + domain = tmp_path / "d.wkt" + domain.write_text(Polygon(ring).wkt) + vtk = run(tmp_path, *mesh_args(quad_dir, domain, "--domain-crs", "EPSG:3035")) + assert_ring_in_output(vtk, to_crs("EPSG:3035", "EPSG:25833", ring)) + assert the_field(vtk, "domain_crs") == "EPSG:3035" + assert the_field(vtk, "domain_transform") == description("EPSG:3035") + + def test_a_crs_with_no_epsg_code(self, tmp_path: Path, quad_dir: Path) -> None: + """16 refused it ("has no EPSG code"); R9 makes `crs` a string for it.""" + ring = to_crs("EPSG:25833", LCC, ACROSS) + domain = tmp_path / "d.wkt" + domain.write_text(Polygon(ring).wkt) + vtk = run(tmp_path, *mesh_args(quad_dir, domain, "--domain-crs", LCC)) + assert_ring_in_output(vtk, to_crs(LCC, "EPSG:25833", ring)) + recorded = the_field(vtk, "domain_crs") + assert recorded.isascii() + assert CRS.from_user_input(recorded) == CRS.from_user_input(LCC) + assert the_field(vtk, "domain_transform") == description(LCC) + + def test_ogc_crs84_member(self, tmp_path: Path, quad_dir: Path) -> None: + lon_lat = to_crs("EPSG:25833", "OGC:CRS84", ACROSS) + domain = write_geojson(tmp_path / "d.geojson", lon_lat, crs="urn:ogc:def:crs:OGC:1.3:CRS84") + vtk = run(tmp_path, *mesh_args(quad_dir, domain)) + assert_ring_in_output(vtk, to_crs("OGC:CRS84", "EPSG:25833", lon_lat)) + assert CRS.from_user_input(the_field(vtk, "domain_crs")) == CRS.from_user_input("OGC:CRS84") + + def test_a_single_file_records_no_tile_list(self, tmp_path: Path, quad_dir: Path) -> None: + inside: Ring = [ + (X0 + 5.3, Y0 - 12.9), + (X0 + 33.7, Y0 - 12.1), + (X0 + 32.9, Y0 - 3.1), + (X0 + 6.1, Y0 - 3.7), + ] + lon_lat = to_crs("EPSG:25833", "EPSG:4326", inside) + domain = write_geojson(tmp_path / "d.geojson", lon_lat) + vtk = run(tmp_path, *mesh_args(quad_dir / "nw.tif", domain)) + assert "dem_tiles" not in vtk.field_data + assert the_field(vtk, "domain_crs") == "EPSG:4326" + + def test_a_domain_that_avoids_a_missing_tile_meshes(self, tmp_path: Path) -> None: + """R4 point 5: NaN filler outside the needed region is never read.""" + source = whole(12, 12, dy=10.0, array=terrain(12, 12)) + write_tiles(tmp_path / "ell", blocks(source, 6, 6, skip=[(1, 1)])) + ell: Ring = [ + (X0 + 5.5, Y0 - 5.5), + (X0 + 105.5, Y0 - 5.5), + (X0 + 105.5, Y0 - 35.0), + (X0 + 35.0, Y0 - 35.0), + (X0 + 35.0, Y0 - 105.5), + (X0 + 5.5, Y0 - 105.5), + ] + lon_lat = to_crs("EPSG:25833", "EPSG:4326", ell) + domain = write_geojson(tmp_path / "ell.geojson", lon_lat) + vtk = run(tmp_path, *mesh_args(tmp_path / "ell", domain)) + assert np.isfinite(vtk.points).all() + assert_inside(vtk, to_crs("EPSG:4326", "EPSG:25833", lon_lat)) + + +class TestTheSameCrs: + """ "A domain already in the DEM's CRS is not transformed, and the mesh is + bit-identical to 16's." Relational, on one machine (module docstring): the + run as it is against the run with 16's data flow put back.""" + + @pytest.fixture + def bumpy(self, tmp_path: Path) -> Path: + array = np.random.default_rng(16).uniform(0.0, 50.0, (ROWS, COLS)).astype(np.float32) + return write_tiff(tmp_path / "bumpy.tif", micro_tiff(array)) + + @pytest.fixture + def no_transformer(self, monkeypatch: pytest.MonkeyPatch) -> None: + def refuse(*args: Any, **kwargs: Any) -> Any: + raise AssertionError(f"Transformer.from_crs{args} was called") + + monkeypatch.setattr(Transformer, "from_crs", staticmethod(refuse)) + + @pytest.fixture + def handed_on(self, monkeypatch: pytest.MonkeyPatch) -> list[DomainPolygon | None]: + """Every domain the CLI hands `_dem_mesh`, in call order.""" + seen: list[DomainPolygon | None] = [] + real = cli._dem_mesh + + def spy(*args: Any, **kwargs: Any) -> Any: + seen.append(kwargs.get("domain", args[7] if len(args) > 7 else None)) + return real(*args, **kwargs) + + monkeypatch.setattr(cli, "_dem_mesh", spy) + return seen + + @staticmethod + def as_16_opened_it(request: DemRequest) -> DemInput: + """16's data flow for one file and a domain: the whole file, opened as + without a domain, and the domain object as `read_domain` returned it. + Neither `DomainPolygon.to_crs` nor the window cut to the domain runs.""" + opened = open_dem(request.model_copy(update={"domain": None})) + return dataclasses.replace(opened, domain=request.domain) + + def assert_as_16( + self, + tmp_path: Path, + monkeypatch: pytest.MonkeyPatch, + handed_on: list[DomainPolygon | None], + domain: Path, + *args: str, + ) -> None: + now = run(tmp_path, *args, out="now.vtk") + with monkeypatch.context() as m: + m.setattr(cli, "open_dem", self.as_16_opened_it) + before = run(tmp_path, *args, out="before.vtk") + assert len(handed_on) == 2, handed_on + read = read_domain(domain) + for given in handed_on: + assert given is not None + assert given.crs == read.crs + for got, want in zip(rings(given), rings(read), strict=True): + assert got.dtype == want.dtype and got.tobytes() == want.tobytes() + assert len(now.points) > 0 + assert digest(now) == digest(before) + + def test_the_mesh_is_16s_bit_for_bit( + self, + tmp_path: Path, + bumpy: Path, + no_transformer: None, + handed_on: list[DomainPolygon | None], + monkeypatch: pytest.MonkeyPatch, + ) -> None: + square = utm33_geojson(tmp_path / "square.geojson", SQUARE, (HOLE,)) + self.assert_as_16(tmp_path, monkeypatch, handed_on, square, *mesh_args(bumpy, square)) + + def test_the_fields_say_so(self, tmp_path: Path, bumpy: Path) -> None: + square = utm33_geojson(tmp_path / "square.geojson", SQUARE, (HOLE,)) + vtk = run(tmp_path, *mesh_args(bumpy, square)) + assert the_field(vtk, "domain_crs") == "EPSG:25833" + assert "none" in the_field(vtk, "domain_transform").lower() + + @needs_codecs + def test_the_quarter_circle_is_16s_bit_for_bit( + self, + tmp_path: Path, + no_transformer: None, + handed_on: list[DomainPolygon | None], + monkeypatch: pytest.MonkeyPatch, + ) -> None: + domain = utm33_geojson(tmp_path / "quarter.geojson", quarter_circle()) + args = ("--dem", str(KARTVERKET), "--domain", str(domain), "--tolerance", "10") + self.assert_as_16(tmp_path, monkeypatch, handed_on, domain, *args) + + +def rings(domain: DomainPolygon) -> list[np.ndarray]: + p = domain.polygon + return [np.asarray(r.coords) for r in (p.exterior, *p.interiors)] + + +class TestUsage: + def test_bbox_with_domain_is_refused(self, tmp_path: Path, quad_dir: Path) -> None: + """R11: `--bbox` excludes `--domain` [15b]. The domain is in the DEM's + own CRS, so nothing but the flag pair can refuse it.""" + domain = utm33_geojson(tmp_path / "d.geojson", ACROSS) + refused( + tmp_path, + *mesh_args(quad_dir, domain), + "--bbox", "500000", "6599960", "500120", "6600000", + says=("--bbox", "--domain"), + ) # fmt: skip + + @pytest.mark.parametrize("bad", ["EPSG:999999", "not-a-crs"]) + def test_an_unknown_domain_crs_is_named(self, tmp_path: Path, quad_dir: Path, bad: str) -> None: + domain = tmp_path / "d.wkt" + domain.write_text(Polygon(ACROSS).wkt) + refused(tmp_path, *mesh_args(quad_dir, domain, "--domain-crs", bad), says=(bad,)) + + +class TestRefusedAfterTheTransform: + """R9 and "Degeneracy policy": what lands outside the DEM is refused by the + extent check, naming the domain's CRS and the DEM's.""" + + def test_a_domain_far_from_the_dem(self, tmp_path: Path, quad_dir: Path) -> None: + far = [(x - 200_000.0, y + 50_000.0) for x, y in ACROSS] + domain = write_geojson(tmp_path / "d.geojson", to_crs("EPSG:25833", "EPSG:4326", far)) + refused(tmp_path, *mesh_args(quad_dir, domain), says=("4326", "25833")) + + def test_a_domain_reaching_past_the_tiles(self, tmp_path: Path, quad_dir: Path) -> None: + past = [*ACROSS[:2], (X0 + 123.0, Y0 - 3.7), ACROSS[3]] + domain = write_geojson(tmp_path / "d.geojson", to_crs("EPSG:25833", "EPSG:4326", past)) + refused(tmp_path, *mesh_args(quad_dir, domain), says=("--domain", "outside", "4326")) + + def test_latitude_first_is_a_wrong_file(self, tmp_path: Path, quad_dir: Path) -> None: + """ "A GeoJSON with latitude first is a wrong file, not a case, and is + refused by the extent check." """ + lon_lat = to_crs("EPSG:25833", "EPSG:4326", ACROSS) + domain = write_geojson(tmp_path / "d.geojson", [(y, x) for x, y in lon_lat]) + refused(tmp_path, *mesh_args(quad_dir, domain), says=("4326", "25833")) + + def test_axis_order_is_able_to_fail( + self, tmp_path: Path, quad_dir: Path, monkeypatch: pytest.MonkeyPatch + ) -> None: + """The same run as `test_wgs84_geojson_without_a_crs_member`, with every + transformer built latitude first: the vertices land thousands of km + away, and the extent check must refuse them.""" + real = Transformer.from_crs + + def latitude_first(*args: Any, **kwargs: Any) -> Any: + kwargs["always_xy"] = False + return real(*args, **kwargs) + + lon_lat = to_crs("EPSG:25833", "EPSG:4326", ACROSS) + domain = write_geojson(tmp_path / "catchment.geojson", lon_lat) + run(tmp_path, *mesh_args(quad_dir, domain)) # the control: unpatched, it meshes + monkeypatch.setattr(Transformer, "from_crs", staticmethod(latitude_first)) + refused(tmp_path, *mesh_args(quad_dir, domain), says=("4326", "25833")) + + +class TestSeamsWithADomain: + """Ola's Q1 revised through `--domain`: `dem_seams` is recorded as without + one. `test_cli_mesh_mosaic.TestQ1Seams.disagreeing` (reached through its + module, so pytest does not collect that class twice) plants +4 at global node (2, 6), which is + (500 060, 6 599 990), inside `ACROSS`.""" + + def test_a_disagreeing_pair_is_recorded(self, tmp_path: Path) -> None: + source = whole(9, 13, array=terrain(9, 13)) + dem = test_cli_mesh_mosaic.TestQ1Seams.disagreeing(tmp_path, source) + lon_lat = to_crs("EPSG:25833", "EPSG:4326", ACROSS) + vtk = run(tmp_path, *mesh_args(dem, write_geojson(tmp_path / "c.geojson", lon_lat))) + assert the_field(vtk, "domain_crs") == "EPSG:4326" + assert the_field(vtk, "dem_seams") == "ne.tif | nw.tif: nodes 1, max 4, median 4" + + def test_agreeing_tiles_record_none(self, tmp_path: Path, quad_dir: Path) -> None: + lon_lat = to_crs("EPSG:25833", "EPSG:4326", ACROSS) + vtk = run(tmp_path, *mesh_args(quad_dir, write_geojson(tmp_path / "c.geojson", lon_lat))) + assert the_field(vtk, "dem_seams") == "none" + + @pytest.mark.parametrize( + ("shape", "recorded"), + [("STRIP", "none"), ("ELL", "ne.tif | nw.tif: nodes 2, max 2, median 1.25")], + ) + def test_only_the_needed_region_is_counted( + self, tmp_path: Path, shape: str, recorded: str + ) -> None: + """Ola, 2026-09-28. `test_dem_input_domain.TestSeamsInsideTheNeededRegion` + (reached through its module, so pytest does not collect it twice): the + strip's region misses every planted node, the L's takes two of four.""" + planted = test_dem_input_domain.TestSeamsInsideTheNeededRegion + dem = planted.disagreeing(tmp_path) + lon_lat = to_crs("EPSG:25833", "EPSG:4326", planted.utm33(getattr(planted, shape))) + vtk = run(tmp_path, *mesh_args(dem, write_geojson(tmp_path / "c.geojson", lon_lat))) + assert the_field(vtk, "dem_seams") == recorded + + +class TestRealSeam: + """The 15b Acceptance on the committed DTM10 extract: a catchment-like + polygon in EPSG:4326 over the 6400_4 | 6400_1 seam, across its overlap.""" + + UTM33_RING: ClassVar[Ring] = [ + (47_503.3, 6_467_207.7), + (52_496.1, 6_467_301.9), + (52_402.7, 6_469_298.3), + (47_601.9, 6_469_203.1), + ] + + def test_meshes_with_the_domain_fields(self, tmp_path: Path) -> None: + seam = Path(__file__).resolve().parents[1] / "fixtures" / "dtm10" / "seam" + lon_lat = to_crs("EPSG:25833", "EPSG:4326", self.UTM33_RING) + domain = write_geojson(tmp_path / "catchment.geojson", lon_lat) + vtk = run(tmp_path, *mesh_args(seam, domain)) + assert the_field(vtk, "dem_tiles") == "6400_1_10m_z33.tif; 6400_4_10m_z33.tif" + assert the_field(vtk, "domain_crs") == "EPSG:4326" + assert the_field(vtk, "domain_transform") == description("EPSG:4326") + assert_inside(vtk, to_crs("EPSG:4326", "EPSG:25833", lon_lat)) + assert np.isfinite(vtk.points).all() diff --git a/tests/python/test_crs.py b/tests/python/test_crs.py new file mode 100644 index 00000000..ebf15c12 --- /dev/null +++ b/tests/python/test_crs.py @@ -0,0 +1,196 @@ +"""`tin_engine.crs`: `parse_crs` and `reprojector`, the one transform site (increment 15b). + +`docs/increments/15-dem-mosaic.md` R2 ([15b] `tin_engine/crs.py`), R9 and I8. +Pinned by this suite (see "Pinned by the red suite (15b)"): + +- `parse_crs(text) -> pyproj.CRS` accepts anything `CRS.from_user_input` + accepts, and refuses anything else with a `ValueError` naming the text. + pyproj's own `CRSError` is a `RuntimeError`, not a `ValueError`, so letting it + through would escape every `except ValueError` that turns a refusal into a + usage error. +- `reprojector(src, dst)` takes two CRSs (text or `pyproj.CRS`) and returns a + callable mapping an `(N, 2)` array-like of `(x, y)` to an `(N, 2)` float64 + array. `x` is easting or longitude, `y` northing or latitude, whatever the + CRS's own axis order (`always_xy`), and the result is pyproj's own + `always_xy` transform bit for bit. +- The only `Transformer.from_crs` call in `src_python/` is in + `tin_engine/crs.py` (I8). `test_always_xy.py` still guards that the call + passes a literal `always_xy=True`. + +Real coordinates are checked against pyproj directly, never against numbers +typed into this file. The oracle is `Transformer.from_crs(..., always_xy=True)` +built here, in the test, which the grep guard does not scan. + +HOW THIS FILE GOES RED: `tin_engine.crs` is imported inside a fixture, so each +test fails on its own with `ModuleNotFoundError` and collection is unaffected. +""" + +from __future__ import annotations + +import ast +import importlib +from pathlib import Path +from types import ModuleType + +import numpy as np +import pytest +from pyproj import CRS, Transformer + +SRC_PYTHON = Path(__file__).resolve().parents[2] / "src_python" + +# Points in (longitude, latitude): the committed DTM10 seam extract near +# Flekkefjord, the benchmark tile's corner in Finnmark, Oslo, and the UTM 33 +# central meridian. Longitude and latitude differ in every row, so swapping +# them cannot give the same answer. +LON_LAT = np.array( + [ + [7.35061409316692, 58.12324881275789], + [21.8, 71.1], + [10.75, 59.91], + [15.0, 59.5], + ] +) + + +@pytest.fixture(scope="module") +def crs() -> ModuleType: + return importlib.import_module("tin_engine.crs") + + +def oracle(src: str, dst: str, xy: np.ndarray) -> np.ndarray: + t = Transformer.from_crs(src, dst, always_xy=True) + x, y = t.transform(xy[:, 0], xy[:, 1]) + return np.column_stack([x, y]) + + +class TestParseCrs: + @pytest.mark.parametrize( + ("text", "epsg"), + [ + ("EPSG:25833", 25833), + ("epsg:25833", 25833), + ("urn:ogc:def:crs:EPSG::25833", 25833), + ("EPSG:4326", 4326), + ], + ) + def test_epsg_spellings(self, crs: ModuleType, text: str, epsg: int) -> None: + parsed = crs.parse_crs(text) + assert isinstance(parsed, CRS) + assert parsed == CRS.from_epsg(epsg) + + def test_a_crs_without_an_epsg_code_is_accepted(self, crs: ModuleType) -> None: + """R9: `crs: str`, since a domain may have no EPSG code.""" + text = "+proj=lcc +lat_1=60 +lat_2=65 +lat_0=62 +lon_0=15 +ellps=GRS80 +units=m +no_defs" + parsed = crs.parse_crs(text) + assert parsed == CRS.from_user_input(text) + assert parsed.to_epsg() is None + + def test_wkt2_is_accepted(self, crs: ModuleType) -> None: + wkt = CRS.from_epsg(3035).to_wkt() + assert crs.parse_crs(wkt) == CRS.from_epsg(3035) + + @pytest.mark.parametrize("text", ["EPSG:999999", "not a crs", "urn:ogc:def:crs:EPSG::0"]) + def test_an_unknown_crs_is_a_value_error_naming_it(self, crs: ModuleType, text: str) -> None: + with pytest.raises(ValueError) as info: + crs.parse_crs(text) + assert text in str(info.value), info.value + + def test_empty_text_is_a_value_error(self, crs: ModuleType) -> None: + with pytest.raises(ValueError): + crs.parse_crs("") + + +class TestReprojector: + def test_geographic_to_utm33_is_pyprojs_always_xy_bit_for_bit(self, crs: ModuleType) -> None: + out = crs.reprojector("EPSG:4326", "EPSG:25833")(LON_LAT) + assert isinstance(out, np.ndarray) + assert out.dtype == np.float64 + assert out.shape == LON_LAT.shape + assert out.tobytes() == oracle("EPSG:4326", "EPSG:25833", LON_LAT).tobytes() + + def test_axis_order_is_x_then_y_for_a_latitude_first_crs(self, crs: ModuleType) -> None: + """EPSG:4326's own axis order is latitude first. The input here is + (longitude, latitude), as in a GeoJSON file and a GeoTIFF's model space, + and the central meridian point lands on UTM 33's false easting.""" + out = crs.reprojector("EPSG:4326", "EPSG:25833")(LON_LAT[3:]) + assert out[0, 0] == pytest.approx(500_000.0, abs=1e-6) + assert 6_500_000.0 < out[0, 1] < 6_700_000.0 + # The same numbers read latitude first land thousands of km away: the + # mutant this test exists to kill. + swapped = Transformer.from_crs("EPSG:4326", "EPSG:25833").transform(*LON_LAT[3]) + assert abs(swapped[0] - out[0, 0]) > 1_000_000.0 + + def test_output_is_x_then_y(self, crs: ModuleType) -> None: + """Northing is the larger number everywhere in Norway, so a column swap + on the way out cannot pass.""" + out = crs.reprojector("EPSG:4326", "EPSG:25833")(LON_LAT) + assert (out[:, 1] > 6_000_000.0).all() + assert (out[:, 0] < 1_200_000.0).all() + + @pytest.mark.parametrize( + ("src", "dst"), + [ + ("EPSG:25832", "EPSG:25833"), + ("EPSG:3035", "EPSG:25833"), + ("EPSG:25833", "EPSG:4326"), + ("OGC:CRS84", "EPSG:25833"), + ], + ) + def test_other_pairs_are_pyprojs_always_xy(self, crs: ModuleType, src: str, dst: str) -> None: + points = oracle("EPSG:4326", src, LON_LAT) + out = crs.reprojector(src, dst)(points) + assert out.tobytes() == oracle(src, dst, points).tobytes() + + def test_crs_objects_are_accepted(self, crs: ModuleType) -> None: + out = crs.reprojector(CRS.from_epsg(4326), crs.parse_crs("EPSG:25833"))(LON_LAT) + assert out.tobytes() == oracle("EPSG:4326", "EPSG:25833", LON_LAT).tobytes() + + def test_a_list_of_pairs_is_accepted(self, crs: ModuleType) -> None: + pairs = [(float(x), float(y)) for x, y in LON_LAT] + out = crs.reprojector("EPSG:4326", "EPSG:25833")(pairs) + assert out.tobytes() == oracle("EPSG:4326", "EPSG:25833", LON_LAT).tobytes() + + def test_no_vertex_is_added_or_dropped(self, crs: ModuleType) -> None: + """R9: vertices only, not densified.""" + for n in (1, 2, 7): + assert crs.reprojector("EPSG:4326", "EPSG:25833")(LON_LAT[:1].repeat(n, 0)).shape == ( + n, + 2, + ) + + +# ---------------------------------------------------------------- I8 + + +def from_crs_sites(root: Path) -> list[str]: + """Every call whose callee is named `from_crs`, as `relative/path.py:line`.""" + found: list[str] = [] + for path in sorted(root.rglob("*.py")): + tree = ast.parse(path.read_text(encoding="utf-8"), filename=str(path)) + for node in ast.walk(tree): + if isinstance(node, ast.Call): + callee = node.func + name = callee.attr if isinstance(callee, ast.Attribute) else None + name = callee.id if isinstance(callee, ast.Name) else name + if name == "from_crs": + found.append(f"{path.relative_to(root).as_posix()}:{node.lineno}") + return found + + +def test_exactly_one_from_crs_site_and_it_is_in_crs_py() -> None: + """I8 and R9: `crs.reprojector` holds the only `Transformer.from_crs` in + `src_python/`. `test_always_xy.py` checks its keywords.""" + sites = from_crs_sites(SRC_PYTHON) + assert len(sites) == 1, sites + assert sites[0].startswith("tin_engine/crs.py:"), sites + + +def test_the_site_finder_finds_planted_sites(tmp_path: Path) -> None: + (tmp_path / "a.py").write_text( + "from pyproj import Transformer\n" + "t = Transformer.from_crs('EPSG:4326', 'EPSG:25833', always_xy=True)\n" + "u = from_crs(1, 2)\n", + encoding="utf-8", + ) + (tmp_path / "b.py").write_text("x = 1\n", encoding="utf-8") + assert from_crs_sites(tmp_path) == ["a.py:2", "a.py:3"] diff --git a/tests/python/test_dem_input_domain.py b/tests/python/test_dem_input_domain.py new file mode 100644 index 00000000..eb6d589a --- /dev/null +++ b/tests/python/test_dem_input_domain.py @@ -0,0 +1,853 @@ +"""`open_dem` with a domain in its own CRS (increment 15b). + +`docs/increments/15-dem-mosaic.md` R1 ("15b adds the domain"), R4 point 5 +(the needed region: "the domain polygon grown by one cell"), R6 (the domain's +bounds in the DEM's CRS take `--bbox`'s place, and the two exclude each other) +and R9 (the extent check runs after the transform, against the mosaic's +coverage). Pinned by this suite (see "Pinned by the red suite (15b)"): + +- `DemRequest(sources=, bounds=, nodata=, domain=)`: `domain` is a + `DomainPolygon` as read, in its own CRS, default `None`. A request with both + `bounds` and `domain` is refused with a `ValueError`. +- `open_dem` moves the domain into the DEM's CRS, plans on its bounds with its + needed region, checks its extent against the plan, and only then assembles. + `DemInput.domain` is the domain in the DEM's CRS (`None` without one). +- The plan is the plan of `bounds` equal to the moved domain's bounds, except + where a domain vertex lies within the 1e-6-cell snap band past a node line: + there the domain's window is one node line wider on that side. A vertex + exactly on a node line, or 2e-6 cell past one, widens nothing (review S1). +- "Grown by one cell" is pinned only away from its edge: a missing node 0.73 + cell from a domain vertex refuses the request; missing nodes 2.5 cells or + more from it are NaN filler. Exactly one cell is not ruled. The cell is the + chosen plan's, so a tile the domain does not select, and the order names + sort in, change nothing (review B1). +- Tiles in more than one CRS are refused with a domain, naming the codes + (review S2). +- Every extent refusal fires before any tile is loaded (I6). +- The seam report (Ola's Q1 revised) survives the domain path: `seams` names + each disagreeing pair with a domain in the DEM's CRS or in another one, and + is `()` when the overlaps agree (test amendment after the 15b review). +- With a domain the seam report counts only nodes inside the needed region, + the domain grown by the plan's cell diagonal, mitred (Ola, 2026-09-28): a + disagreement inside the plan's rectangle but outside that region is not + reported, one just outside the polygon but inside the region is, and which + tile's value a node takes does not change (test amendment for the ruling). + +Synthetic tiles are micro-TIFFs written by `test_dem_input.py`'s helpers; +the real case is the committed DTM10 seam extract (`tests/fixtures/dtm10/`), +with a domain in EPSG:25832, the zone the place lies in. + +HOW THIS FILE GOES RED: `tin_engine.dem_input` and `tin_engine.domain` are +imported inside fixtures, and `DemRequest` has no `domain` field before 15b, +so each test fails on its own and collection is unaffected. +""" + +from __future__ import annotations + +import importlib +import json +import math +from pathlib import Path +from types import ModuleType +from typing import Any, ClassVar + +import numpy as np +import pytest +import shapely +from pyproj import CRS, Transformer +from shapely.geometry import Polygon, box + +from geotiff_fixtures import BASE_KEYS, PROJECTED_CS_TYPE +from mosaic_fixtures import ( + X0, + Y0, + blocks, + deepest_interior, + piece, + quadrants, + same_array, + seams, + seams_of, + whole, +) +from test_dem_input import SEAM, decoded, tiff_of, write_tiles + +Ring = list[tuple[float, float]] +UTM33 = "urn:ogc:def:crs:EPSG::25833" + + +@pytest.fixture(scope="module") +def di() -> ModuleType: + return importlib.import_module("tin_engine.dem_input") + + +@pytest.fixture(scope="module") +def dm() -> ModuleType: + return importlib.import_module("tin_engine.domain") + + +@pytest.fixture(scope="module") +def mz() -> ModuleType: + return importlib.import_module("tin_engine.mosaic") + + +def to_crs(src: str, dst: str, ring: Ring) -> Ring: + """pyproj's own `always_xy` transform of `ring`: the oracle.""" + t = Transformer.from_crs(src, dst, always_xy=True) + xs, ys = t.transform(np.array([p[0] for p in ring]), np.array([p[1] for p in ring])) + return [(float(x), float(y)) for x, y in zip(xs, ys, strict=True)] + + +def domain_file( + path: Path, outer: Ring, holes: tuple[Ring, ...] = (), crs: str | None = UTM33 +) -> Path: + doc: dict[str, Any] = { + "type": "Polygon", + "coordinates": [[*r, r[0]] for r in (outer, *holes)], + } + if crs is not None: + doc["crs"] = {"type": "name", "properties": {"name": crs}} + path.write_text(json.dumps(doc)) + return path + + +def read( + dm: ModuleType, tmp_path: Path, utm33: Ring, crs: str, holes: tuple[Ring, ...] = () +) -> Any: + """A domain file holding the UTM 33 ring `utm33`, written in `crs`. + + `crs` is "EPSG:25833" (written with its `crs` member), "EPSG:4326" (no + member, as RFC 7946 has it) or another EPSG code (with its member).""" + if crs == "EPSG:25833": + return dm.read_domain(domain_file(tmp_path / "d.geojson", utm33, holes)) + outer = to_crs("EPSG:25833", crs, utm33) + moved = tuple(to_crs("EPSG:25833", crs, h) for h in holes) + member = None if crs == "EPSG:4326" else crs + return dm.read_domain(domain_file(tmp_path / "d.geojson", outer, moved, member)) + + +def request(di: ModuleType, *sources: Path, domain: Any = None, bounds: Any = None) -> Any: + return di.DemRequest(sources=tuple(sources), bounds=bounds, domain=domain) + + +@pytest.fixture +def no_load(monkeypatch: pytest.MonkeyPatch) -> None: + """I6: any tile load from here on fails the test.""" + repository = importlib.import_module("tin_engine.io.repository") + + def refuse(self: Any, name: str, *args: Any, **kwargs: Any) -> Any: + raise AssertionError(f"load({name!r}) was called") + + monkeypatch.setattr(repository.TiffDemRepository, "load", refuse) + + +# ------------------------------------------------------------------ tiles + +# `whole(9, 13)` in quadrants, point-registered, one shared node line: +# nodes x 500 000 .. 500 120 (dx 10), y 6 600 000 .. 6 599 960 (dy 5); the +# north-west tile holds rows 0-4 and columns 0-6. +IN_NW: Ring = [ + (X0 + 5.3, Y0 - 12.9), + (X0 + 33.7, Y0 - 12.1), + (X0 + 32.9, Y0 - 3.1), + (X0 + 6.1, Y0 - 3.7), +] +ACROSS_ALL: Ring = [ + (X0 + 12.3, Y0 - 36.3), + (X0 + 107.7, Y0 - 35.9), + (X0 + 106.1, Y0 - 3.7), + (X0 + 13.9, Y0 - 4.1), +] + + +@pytest.fixture +def quad_dir(tmp_path: Path) -> Path: + write_tiles(tmp_path / "quad", quadrants(whole(9, 13), row_cut=4, col_cut=6, overlap=1)) + return tmp_path / "quad" + + +# 12 x 12 nodes, dx = dy = 10, cut into four abutting 6 x 6 blocks with the +# south-east one missing: nodes x 500 060 .. 500 110, y 6 599 940 .. 6 599 890 +# are in no tile. +@pytest.fixture +def l_dir(tmp_path: Path) -> Path: + write_tiles(tmp_path / "ell", blocks(whole(12, 12, dy=10.0), 6, 6, skip=[(1, 1)])) + return tmp_path / "ell" + + +# A triangle whose corner (53.3, -57.1) has the missing node (60, -60) as a +# bilinear corner: 6.7 m east and 2.9 m south, well within one cell. +NEAR_THE_HOLE: Ring = [(X0 + 5.5, Y0 - 5.5), (X0 + 53.3, Y0 - 57.1), (X0 + 5.5, Y0 - 57.1)] +# An L whose box holds every missing node, and whose nearest point to one, +# (35, -35), is 25 m from each axis of (60, -60): 2.5 cells. +AROUND_THE_HOLE: Ring = [ + (X0 + 5.5, Y0 - 5.5), + (X0 + 105.5, Y0 - 5.5), + (X0 + 105.5, Y0 - 35.0), + (X0 + 35.0, Y0 - 35.0), + (X0 + 35.0, Y0 - 105.5), + (X0 + 5.5, Y0 - 105.5), +] + + +class TestDemRequest: + def test_domain_defaults_to_none(self, di: ModuleType, quad_dir: Path) -> None: + assert di.DemRequest(sources=(quad_dir,)).domain is None + + def test_bounds_and_domain_exclude_each_other( + self, di: ModuleType, dm: ModuleType, mz: ModuleType, quad_dir: Path, tmp_path: Path + ) -> None: + """R6 and R11: with a domain, its bounds take `--bbox`'s place.""" + domain = read(dm, tmp_path, IN_NW, "EPSG:25833") + bounds = mz.Bounds(x_min=X0, y_min=Y0 - 20, x_max=X0 + 50, y_max=Y0) + with pytest.raises(ValueError): + request(di, quad_dir, domain=domain, bounds=bounds) + + def test_without_a_domain_the_input_has_none(self, di: ModuleType, quad_dir: Path) -> None: + assert di.open_dem(request(di, quad_dir)).domain is None + + +class TestTheDomainChoosesTheTiles: + @pytest.mark.parametrize("crs", ["EPSG:25833", "EPSG:4326", "EPSG:25832"]) + def test_the_domain_arrives_in_the_dems_crs( + self, di: ModuleType, dm: ModuleType, quad_dir: Path, tmp_path: Path, crs: str + ) -> None: + domain = read(dm, tmp_path, ACROSS_ALL, crs) + opened = di.open_dem(request(di, quad_dir, domain=domain)) + assert CRS.from_user_input(opened.domain.crs) == CRS.from_epsg(25833) + given = list(domain.polygon.exterior.coords) + expected = given if crs == "EPSG:25833" else to_crs(crs, "EPSG:25833", given) + got = {(float(x), float(y)) for x, y in opened.domain.polygon.exterior.coords} + assert got == set(expected) + + @pytest.mark.parametrize("crs", ["EPSG:25833", "EPSG:4326"]) + def test_the_plan_is_the_plan_of_the_moved_domains_bounds( + self, + di: ModuleType, + dm: ModuleType, + mz: ModuleType, + quad_dir: Path, + tmp_path: Path, + crs: str, + ) -> None: + opened = di.open_dem(request(di, quad_dir, domain=read(dm, tmp_path, IN_NW, crs))) + x_min, y_min, x_max, y_max = opened.domain.polygon.bounds + bounds = mz.Bounds(x_min=x_min, y_min=y_min, x_max=x_max, y_max=y_max) + by_box = di.open_dem(request(di, quad_dir, bounds=bounds)) + assert opened.plan == by_box.plan + assert [t.name for t in opened.plan.tiles] == ["nw.tif"] + + def test_a_domain_across_every_tile_selects_every_tile( + self, di: ModuleType, dm: ModuleType, quad_dir: Path, tmp_path: Path + ) -> None: + domain = read(dm, tmp_path, ACROSS_ALL, "EPSG:4326") + opened = di.open_dem(request(di, quad_dir, domain=domain)) + assert [t.name for t in opened.plan.tiles] == ["ne.tif", "nw.tif", "se.tif", "sw.tif"] + + def test_a_single_file_is_cut_to_the_domains_window( + self, di: ModuleType, dm: ModuleType, quad_dir: Path, tmp_path: Path + ) -> None: + """R6 on one file: a window of it, as `--bbox` on one file gives.""" + opened = di.open_dem( + request(di, quad_dir / "nw.tif", domain=read(dm, tmp_path, IN_NW, "EPSG:4326")) + ) + # IN_NW spans columns 0..4 and rows 0..3 once snapped outward. + assert (opened.tile.meta.rows, opened.tile.meta.cols) == (4, 5) + assert opened.tile.meta.x_min == X0 + assert opened.tile.meta.y_max == Y0 + + +class TestNeededRegion: + """R4 point 5 with a domain: the polygon grown by one cell, in the DEM's CRS.""" + + @pytest.mark.parametrize("crs", ["EPSG:25833", "EPSG:4326"]) + def test_a_missing_node_the_boundary_needs_is_refused( + self, + di: ModuleType, + dm: ModuleType, + mz: ModuleType, + l_dir: Path, + tmp_path: Path, + no_load: None, + crs: str, + ) -> None: + domain = read(dm, tmp_path, NEAR_THE_HOLE, crs) + with pytest.raises(mz.MosaicError) as info: + di.open_dem(request(di, l_dir, domain=domain)) + assert "in no tile" in str(info.value) + assert "500060" in str(info.value) + + @pytest.mark.parametrize("crs", ["EPSG:25833", "EPSG:4326"]) + def test_missing_nodes_away_from_the_domain_are_filler( + self, di: ModuleType, dm: ModuleType, l_dir: Path, tmp_path: Path, crs: str + ) -> None: + domain = read(dm, tmp_path, AROUND_THE_HOLE, crs) + opened = di.open_dem(request(di, l_dir, domain=domain)) + tile = opened.tile + assert (tile.meta.rows, tile.meta.cols) == (12, 12) + assert np.isnan(tile.array[6:, 6:]).all() + assert np.isfinite(tile.array[:6, :]).all() + assert np.isfinite(tile.array[:, :6]).all() + + def test_a_domain_enclosing_a_missing_tile_is_refused( + self, di: ModuleType, dm: ModuleType, mz: ModuleType, tmp_path: Path, no_load: None + ) -> None: + """The interior is needed, not only the boundary.""" + write_tiles(tmp_path / "ring", blocks(whole(18, 18, dy=10.0), 6, 6, skip=[(1, 1)])) + around: Ring = [ + (X0 + 5.5, Y0 - 5.5), + (X0 + 164.5, Y0 - 5.5), + (X0 + 164.5, Y0 - 164.5), + (X0 + 5.5, Y0 - 164.5), + ] + domain = read(dm, tmp_path, around, "EPSG:4326") + with pytest.raises(mz.MosaicError) as info: + di.open_dem(request(di, tmp_path / "ring", domain=domain)) + assert "in no tile" in str(info.value) + + def test_a_domain_hole_over_a_missing_tile_is_not_needed( + self, di: ModuleType, dm: ModuleType, tmp_path: Path + ) -> None: + """The domain's own hole, 2.5 cells wider than the missing block on + every side, is not part of the needed region.""" + write_tiles(tmp_path / "ring", blocks(whole(18, 18, dy=10.0), 6, 6, skip=[(1, 1)])) + around: Ring = [ + (X0 + 5.5, Y0 - 5.5), + (X0 + 164.5, Y0 - 5.5), + (X0 + 164.5, Y0 - 164.5), + (X0 + 5.5, Y0 - 164.5), + ] + # Missing nodes: x 60..110, y -60..-110. The hole: x 35..135, y -35..-135. + hole: Ring = [ + (X0 + 35.0, Y0 - 35.0), + (X0 + 35.0, Y0 - 135.0), + (X0 + 135.0, Y0 - 135.0), + (X0 + 135.0, Y0 - 35.0), + ] + domain = read(dm, tmp_path, around, "EPSG:4326", holes=(hole,)) + opened = di.open_dem(request(di, tmp_path / "ring", domain=domain)) + assert np.isnan(opened.tile.array[6:12, 6:12]).all() + + +class TestExtentAfterTheTransform: + """R9: the extent check runs after the transform, against the tiles, and + fires before any tile is loaded (I6).""" + + def test_a_domain_far_from_every_tile_is_refused( + self, di: ModuleType, dm: ModuleType, quad_dir: Path, tmp_path: Path, no_load: None + ) -> None: + far = [(x - 200_000.0, y + 50_000.0) for x, y in IN_NW] + with pytest.raises(ValueError): + di.open_dem(request(di, quad_dir, domain=read(dm, tmp_path, far, "EPSG:4326"))) + + @pytest.mark.parametrize("crs", ["EPSG:25833", "EPSG:4326"]) + def test_a_domain_reaching_past_the_tiles_is_refused_as_outside( + self, + di: ModuleType, + dm: ModuleType, + quad_dir: Path, + tmp_path: Path, + no_load: None, + crs: str, + ) -> None: + """The last node column is x = 500 120; one vertex is 3 m past it. The + window is clamped to the tiles, so the per-node coverage check alone + would not see it.""" + past = [*ACROSS_ALL[:2], (X0 + 123.0, Y0 - 3.7), ACROSS_ALL[3]] + with pytest.raises(ValueError) as info: + di.open_dem(request(di, quad_dir, domain=read(dm, tmp_path, past, crs))) + assert "outside" in str(info.value) + + def test_a_domain_on_the_last_node_line_is_inside( + self, di: ModuleType, dm: ModuleType, quad_dir: Path, tmp_path: Path + ) -> None: + """16's border rule, unchanged: the node rectangle is closed.""" + corners: Ring = [(X0, Y0 - 40.0), (X0 + 120.0, Y0 - 40.0), (X0 + 120.0, Y0), (X0, Y0)] + opened = di.open_dem( + request(di, quad_dir, domain=read(dm, tmp_path, corners, "EPSG:25833")) + ) + assert (opened.tile.meta.rows, opened.tile.meta.cols) == (9, 13) + + +class TestRealSeam: + """The committed DTM10 seam (6400_4 | 6400_1, EPSG:25833, near Flekkefjord, + which lies in UTM zone 32), with a domain drawn in EPSG:25832.""" + + # Across the 51-column overlap at x 49 750 .. 50 250; every vertex off-node. + ACROSS_THE_SEAM: ClassVar[Ring] = [ + (47_503.3, 6_467_207.7), + (52_496.1, 6_467_301.9), + (52_402.7, 6_469_298.3), + (47_601.9, 6_469_203.1), + ] + + def test_the_extract_is_where_this_test_thinks(self) -> None: + node_rect = box(47_190.0, 6_466_980.0, 52_810.0, 6_469_530.0) + assert node_rect.contains(Polygon(self.ACROSS_THE_SEAM)) + + def test_a_utm32_domain_opens_both_tiles( + self, di: ModuleType, dm: ModuleType, tmp_path: Path + ) -> None: + domain = read(dm, tmp_path, self.ACROSS_THE_SEAM, "EPSG:25832") + opened = di.open_dem(request(di, SEAM, domain=domain)) + assert [t.name for t in opened.plan.tiles] == ["6400_1_10m_z33.tif", "6400_4_10m_z33.tif"] + given = list(domain.polygon.exterior.coords) + got = {(float(x), float(y)) for x, y in opened.domain.polygon.exterior.coords} + assert got == set(to_crs("EPSG:25832", "EPSG:25833", given)) + + +class TestNeededRegionIsGrownByThePlansSpacing: + """R4 point 5: "grown by one cell" is one cell of the lattice the plan is + on, so the outcome does not depend on a tile the domain does not select, + nor on where that tile's name sorts (15b review, B1). Each repository holds + two spacings in one EPSG: an L of blocks with its south-east block missing, + and one far-away tile on the other spacing, named to sort first or last.""" + + FAR = 5_000.0 # the far tile's west edge, metres east of X0 + + @staticmethod + def far_fine() -> Any: + """A 1 m tile 5 km east of every coarse node: no domain here meets it.""" + return whole(4, 4, dx=1.0, dy=1.0, x_min=X0 + TestNeededRegionIsGrownByThePlansSpacing.FAR) + + @staticmethod + def far_coarse() -> Any: + return whole( + 4, 4, dx=10.0, dy=10.0, x_min=X0 + TestNeededRegionIsGrownByThePlansSpacing.FAR + ) + + @pytest.mark.parametrize("fine", [None, "a_fine.tif", "z_fine.tif"]) + def test_a_node_one_coarse_cell_away_is_needed_whatever_else_is_listed( + self, + di: ModuleType, + dm: ModuleType, + mz: ModuleType, + tmp_path: Path, + no_load: None, + fine: str | None, + ) -> None: + """`NEAR_THE_HOLE` on the 10 m L: the missing node (60, -60) is 7.3 m + from the domain, inside one 10 m cell and outside one 1 m cell. The + plan is on the 10 m lattice, so the request is refused, with or + without a 1 m tile listed first or last.""" + tiles = blocks(whole(12, 12, dy=10.0), 6, 6, skip=[(1, 1)]) + if fine is not None: + tiles[fine] = self.far_fine() + write_tiles(tmp_path / "dem", tiles) + domain = read(dm, tmp_path, NEAR_THE_HOLE, "EPSG:25833") + with pytest.raises(mz.MosaicError) as info: + di.open_dem(request(di, tmp_path / "dem", domain=domain)) + assert "1 nodes the request needs are in no tile" in str(info.value) + assert "500060" in str(info.value) + + @pytest.mark.parametrize("coarse", [None, "a_coarse.tif", "z_coarse.tif"]) + def test_a_node_past_one_fine_cell_is_not_needed_whatever_else_is_listed( + self, di: ModuleType, dm: ModuleType, tmp_path: Path, coarse: str | None + ) -> None: + """The mirror, on a 1 m L (missing nodes x 6..11, y -6..-11) with a far + 10 m tile: `AROUND_THE_HOLE` scaled by a tenth is 2.5 m per axis + (3.5 m) from the nearest missing node, past one 1 m cell and inside one + 10 m cell. The plan is on the 1 m lattice, so the missing block is NaN + filler, with or without the 10 m tile listed first or last.""" + tiles = blocks(whole(12, 12, dx=1.0, dy=1.0), 6, 6, skip=[(1, 1)]) + if coarse is not None: + tiles[coarse] = self.far_coarse() + write_tiles(tmp_path / "dem", tiles) + tenth = [(X0 + (x - X0) / 10, Y0 + (y - Y0) / 10) for x, y in AROUND_THE_HOLE] + opened = di.open_dem( + request(di, tmp_path / "dem", domain=read(dm, tmp_path, tenth, "EPSG:25833")) + ) + assert opened.plan.meta.delta_x == 1.0 + assert (opened.tile.meta.rows, opened.tile.meta.cols) == (12, 12) + assert np.isnan(opened.tile.array[6:, 6:]).all() + assert np.isfinite(opened.tile.array[:6, :]).all() + assert np.isfinite(opened.tile.array[:, :6]).all() + + @pytest.mark.parametrize("fine", ["a_fine.tif", "z_fine.tif"]) + def test_an_unselected_tile_does_not_change_the_plan( + self, di: ModuleType, dm: ModuleType, tmp_path: Path, fine: str + ) -> None: + write_tiles(tmp_path / "alone", blocks(whole(12, 12, dy=10.0), 6, 6, skip=[(1, 1)])) + tiles = blocks(whole(12, 12, dy=10.0), 6, 6, skip=[(1, 1)]) + tiles[fine] = self.far_fine() + write_tiles(tmp_path / "with", tiles) + domain = read(dm, tmp_path, AROUND_THE_HOLE, "EPSG:25833") + alone = di.open_dem(request(di, tmp_path / "alone", domain=domain)) + with_fine = di.open_dem(request(di, tmp_path / "with", domain=domain)) + assert with_fine.plan == alone.plan + + # A rectangle x 5.5..53.3, y -5.5..-57.1 with a notch x 20..30 down to + # y -20 cut from its north edge. Its south-east corner is `NEAR_THE_HOLE`'s, + # 7.3 m from the missing 10 m node (60, -60). + NOTCHED: ClassVar[Ring] = [ + (X0 + 5.5, Y0 - 5.5), + (X0 + 20.0, Y0 - 5.5), + (X0 + 20.0, Y0 - 20.0), + (X0 + 30.0, Y0 - 20.0), + (X0 + 30.0, Y0 - 5.5), + (X0 + 53.3, Y0 - 5.5), + (X0 + 53.3, Y0 - 57.1), + (X0 + 5.5, Y0 - 57.1), + ] + + @staticmethod + def fine(rows: int, cols: int, west: float, north: float) -> Any: + return whole(rows, cols, dx=1.0, dy=1.0, x_min=X0 + west, y_max=Y0 + north) + + @pytest.mark.parametrize("prefix", ["a_", "z_"]) + def test_the_needed_region_grows_again_when_the_lattice_changes( + self, di: ModuleType, dm: ModuleType, mz: ModuleType, tmp_path: Path, prefix: str + ) -> None: + """The re-plan loop runs until the reach is the chosen lattice's cell. + + The 10 m L plus four 1 m tiles covering `NOTCHED` except a 9 x 15-node + gap at x 21..29, y -5..-19, which lies 1 m outside the notch. On the + polygon itself the 1 m lattice is chosen (four tiles against three); + grown by its 1.41 m diagonal the gap is needed, so the 10 m lattice is + chosen; grown by that one's 14.1 m diagonal, (60, -60) is needed, so + the 10 m lattice no longer covers either, and the request is refused + as mixed-lattice (Q5). The 1 m tiles are selected at every stage: + selection follows the ungrown box. Growing only + once, by the first plan's cell, accepts it with (60, -60) as NaN + filler: the B1 bug by another route.""" + tiles = blocks(whole(12, 12, dy=10.0), 6, 6, skip=[(1, 1)]) + tiles[f"{prefix}left.tif"] = self.fine(54, 16, 5.0, -5.0) + tiles[f"{prefix}mid.tif"] = self.fine(39, 9, 21.0, -20.0) + tiles[f"{prefix}rtop.tif"] = self.fine(26, 25, 30.0, -5.0) + tiles[f"{prefix}rbot.tif"] = self.fine(28, 25, 30.0, -31.0) + write_tiles(tmp_path / "dem", tiles) + domain = read(dm, tmp_path, self.NOTCHED, "EPSG:25833") + with pytest.raises(mz.MosaicError) as info: + di.open_dem(request(di, tmp_path / "dem", domain=domain)) + assert "the request selects tiles on two lattices" in str(info.value) + assert f"{prefix}left.tif" in str(info.value) + + +class TestTheSnapBand: + """The domain's window against `--bbox`'s, on `quad_dir` (node lines every + 10 m in x, 5 m in y). 15a's window snaps a box edge within 1e-6 cell past + a node line onto it; a domain vertex there still needs the next line, so + the domain's window is one line wider on that side (15b, `_past`). A + vertex exactly on a node line, or 2e-6 cell past one, adds nothing (15b + review, S1: `<` made `<=` would add a line for every vertex on a node + line, which a gridded catchment has on every side).""" + + # Every edge on a node line, interior to the 9 x 13 node rectangle: + # columns 2..6, rows 2..6. + ON_LINES = (X0 + 20.0, Y0 - 30.0, X0 + 60.0, Y0 - 10.0) + # Unit vector out of the rectangle for each of x_min, y_min, x_max, y_max, + # and the (rows, cols) and (x_min, y_max) change one more line there makes. + OUT = (-1, -1, 1, 1) + + def plans( + self, + di: ModuleType, + dm: ModuleType, + mz: ModuleType, + quad_dir: Path, + tmp_path: Path, + edges: tuple[float, float, float, float], + ) -> tuple[Any, Any]: + """The domain's plan and `--bbox`'s, for the rectangle `edges`.""" + x0, y0, x1, y1 = edges + ring: Ring = [(x0, y0), (x1, y0), (x1, y1), (x0, y1)] + domain = read(dm, tmp_path, ring, "EPSG:25833") + assert domain.polygon.bounds == edges # bit for bit: same CRS, no transform + by_domain = di.open_dem(request(di, quad_dir, domain=domain)).plan + bounds = mz.Bounds(x_min=x0, y_min=y0, x_max=x1, y_max=y1) + by_box = di.open_dem(request(di, quad_dir, bounds=bounds)).plan + return by_domain, by_box + + def test_on_node_lines_the_domain_window_is_the_bbox_window( + self, di: ModuleType, dm: ModuleType, mz: ModuleType, quad_dir: Path, tmp_path: Path + ) -> None: + by_domain, by_box = self.plans(di, dm, mz, quad_dir, tmp_path, self.ON_LINES) + assert (by_box.meta.rows, by_box.meta.cols) == (5, 5) + assert (by_box.meta.x_min, by_box.meta.y_max) == (X0 + 20.0, Y0 - 10.0) + assert by_domain == by_box + + @pytest.mark.parametrize("edge", [0, 1, 2, 3], ids=["x_min", "y_min", "x_max", "y_max"]) + @pytest.mark.parametrize("cells", [0.0, 2e-6]) + def test_outside_the_band_the_domain_window_is_the_bbox_window( + self, + di: ModuleType, + dm: ModuleType, + mz: ModuleType, + quad_dir: Path, + tmp_path: Path, + edge: int, + cells: float, + ) -> None: + edges = list(self.ON_LINES) + edges[edge] += self.OUT[edge] * cells * (10.0 if edge % 2 == 0 else 5.0) + by_domain, by_box = self.plans(di, dm, mz, quad_dir, tmp_path, tuple(edges)) + assert by_domain == by_box + + @pytest.mark.parametrize("edge", [0, 1, 2, 3], ids=["x_min", "y_min", "x_max", "y_max"]) + def test_in_the_band_the_domain_window_is_one_line_wider( + self, + di: ModuleType, + dm: ModuleType, + mz: ModuleType, + quad_dir: Path, + tmp_path: Path, + edge: int, + ) -> None: + """A vertex 1e-7 cell past a node line: `--bbox` snaps onto the line, + the domain's window takes the next one on that side only.""" + edges = list(self.ON_LINES) + step = 10.0 if edge % 2 == 0 else 5.0 + edges[edge] += self.OUT[edge] * 1e-7 * step + by_domain, by_box = self.plans(di, dm, mz, quad_dir, tmp_path, tuple(edges)) + exact, _ = self.plans(di, dm, mz, quad_dir, tmp_path, self.ON_LINES) + assert by_box == exact + m, e = by_domain.meta, exact.meta + wider = { + 0: (e.rows, e.cols + 1, e.x_min - 10.0, e.y_max), + 1: (e.rows + 1, e.cols, e.x_min, e.y_max), + 2: (e.rows, e.cols + 1, e.x_min, e.y_max), + 3: (e.rows + 1, e.cols, e.x_min, e.y_max + 5.0), + }[edge] + assert (m.rows, m.cols, m.x_min, m.y_max) == wider + + +class TestOneCrs: + def test_tiles_in_two_crss_with_a_domain_are_refused_before_any_load( + self, di: ModuleType, dm: ModuleType, mz: ModuleType, tmp_path: Path, no_load: None + ) -> None: + """R6: the domain is moved into the DEM's CRS, so the DEM must have one + (15b review, S2). Both tiles are in the domain's box.""" + tiles = quadrants(whole(9, 13), row_cut=4, col_cut=6, overlap=1) + write_tiles(tmp_path / "two", {"nw.tif": tiles["nw.tif"]}) + utm32 = {**BASE_KEYS, PROJECTED_CS_TYPE: 25832} + (tmp_path / "two" / "ne.tif").write_bytes( + tiff_of(tiles["ne.tif"], geokeys=utm32).getvalue() + ) + domain = read(dm, tmp_path, IN_NW, "EPSG:25833") + with pytest.raises(mz.MosaicError) as info: + di.open_dem(request(di, tmp_path / "two", domain=domain)) + message = str(info.value) + assert "the tiles are in 2 CRSs" in message + assert "EPSG:[25832, 25833]" in message + assert "a domain needs one" in message + + +class TestSeamsOnTheDomainPath: + """Ola's Q1 revised, with a domain: `DemInput.seams` is the mosaic's report, + as without one. `ne.tif` is planted off `nw.tif` on their shared column + (global column 6, in both tiles only): +0.5 at row 1, +2.0 at row 2, both + at least 1 mm, and +0.0005 at row 3, below it. So `ne.tif | nw.tif`, + nodes 2, max 2.0, median 1.25. Every planted node, (500 060, y 6 599 995 + .. 6 599 985), lies inside `ACROSS_ALL`. A resolution that drops the seams + whenever a domain is given fails the disagreeing cases (15b review).""" + + @staticmethod + def disagreeing(tmp_path: Path) -> Path: + source = whole(9, 13) + tiles = quadrants(source, row_cut=4, col_cut=6, overlap=1) + changed = np.array(tiles["ne.tif"].array) + for row, by in ((1, 0.5), (2, 2.0), (3, 0.0005)): + changed[row, 0] += np.float32(by) + tiles["ne.tif"] = piece(source, 0, 5, 6, 13, array=changed) + write_tiles(tmp_path / "disagree", tiles) + return tmp_path / "disagree" + + @pytest.mark.parametrize("crs", ["EPSG:25833", "EPSG:4326", "EPSG:25832"]) + def test_a_disagreeing_pair_is_reported_with_a_domain( + self, di: ModuleType, dm: ModuleType, tmp_path: Path, crs: str + ) -> None: + dem = self.disagreeing(tmp_path) + domain = read(dm, tmp_path, ACROSS_ALL, crs) + opened = di.open_dem(request(di, dem, domain=domain)) + assert opened.domain is not None + assert seams(opened) == [("ne.tif", "nw.tif", 2, 2.0, 1.25)] + + def test_the_report_is_the_one_without_a_domain( + self, di: ModuleType, dm: ModuleType, mz: ModuleType, tmp_path: Path + ) -> None: + """The same report as the node-by-node oracle over the domain's mosaic, + and as `--bbox` at the moved domain's bounds (the same plan, R6).""" + dem = self.disagreeing(tmp_path) + domain = read(dm, tmp_path, ACROSS_ALL, "EPSG:4326") + opened = di.open_dem(request(di, dem, domain=domain)) + tiles = {p.name: decoded(p) for p in sorted(dem.iterdir())} + assert seams(opened) == seams_of(tiles, opened.tile.meta) + x_min, y_min, x_max, y_max = opened.domain.polygon.bounds + bounds = mz.Bounds(x_min=x_min, y_min=y_min, x_max=x_max, y_max=y_max) + assert seams(opened) == seams(di.open_dem(request(di, dem, bounds=bounds))) + + @pytest.mark.parametrize("crs", ["EPSG:25833", "EPSG:4326"]) + def test_agreeing_overlaps_report_nothing_with_a_domain( + self, di: ModuleType, dm: ModuleType, quad_dir: Path, tmp_path: Path, crs: str + ) -> None: + opened = di.open_dem(request(di, quad_dir, domain=read(dm, tmp_path, ACROSS_ALL, crs))) + assert [t.name for t in opened.plan.tiles] == ["ne.tif", "nw.tif", "se.tif", "sw.tif"] + assert opened.seams == () + + +class TestSeamsInsideTheNeededRegion: + """Ola, 2026-09-28 ("yes, go with a"): with a domain, the seam report counts + only nodes inside the needed region, the domain grown by the chosen plan's + cell diagonal, mitred (`_domain_plan`). Without one it is unchanged, and + which tile's value a node takes is unchanged everywhere. + + `whole(21, 21, dy=10)` in quadrants, one shared node line at each cut: + nodes x 500 000 .. 500 200, y 6 600 000 .. 6 599 800, dx = dy = 10, so the + cell diagonal is 14.14 m. `ne.tif` is planted off `nw.tif` on their shared + column x 500 100 (in both tiles only) at y 6 599 990 .. 6 599 950: + +0.5, +2.0, +0.0005 (below 1 mm), +1.0, +4.0. Over all of them the pair is + nodes 4, max 4.0, median 1.5; over the two northmost, nodes 2, max 2.0, + median 1.25. Local coordinates below are metres east of X0 and north of + Y0. Every vertex sits off the node lines and off the snap band, so each + domain's plan is `--bbox`'s at its bounds (checked, not assumed).""" + + PLANTED: ClassVar[tuple[tuple[int, float], ...]] = ( + (1, 0.5), + (2, 2.0), + (3, 0.0005), + (4, 1.0), + (5, 4.0), + ) + ALL: ClassVar[list[tuple[str, str, int, float, float]]] = [("ne.tif", "nw.tif", 4, 4.0, 1.5)] + NORTH_TWO: ClassVar[list[tuple[str, str, int, float, float]]] = [ + ("ne.tif", "nw.tif", 2, 2.0, 1.25) + ] + # A thin strip from the south-west corner to the north-east one. It + # crosses the shared column at y -100 (where the tiles agree); the nearest + # planted node, (100, -50), is 50 / sqrt 2 = 35.4 m from its centre line + # and 31.9 m from its edge, more than twice the 14.14 m growth. + STRIP: ClassVar[list[tuple[float, float]]] = [ + (5.5, -195.5), + (10.5, -195.5), + (195.5, -10.5), + (195.5, -5.5), + (190.5, -5.5), + (5.5, -190.5), + ] + # An L: a north arm y -5.5 .. -12.5 across the whole width, crossing the + # shared column, and a west arm x 5.5 .. 15.5 the whole height. Grown by + # 14.14 m the north arm reaches y -26.6: (100, -10) is inside the polygon, + # (100, -20) outside it and inside the region, (100, -40) and (100, -50) + # outside both, and inside the plan's rectangle. + ELL: ClassVar[list[tuple[float, float]]] = [ + (5.5, -5.5), + (195.5, -5.5), + (195.5, -12.5), + (15.5, -12.5), + (15.5, -195.5), + (5.5, -195.5), + ] + # Wholly west of the shared column: its east edge is 7 m short of it, so + # (100, -10) and (100, -20) are outside the polygon, within one diagonal. + WEST: ClassVar[list[tuple[float, float]]] = [ + (60.5, -5.5), + (93.0, -5.5), + (93.0, -25.5), + (60.5, -25.5), + ] + + @classmethod + def disagreeing(cls, tmp_path: Path) -> Path: + source = whole(21, 21, dy=10.0) + tiles = quadrants(source, row_cut=10, col_cut=10, overlap=1) + changed = np.array(tiles["ne.tif"].array) + for row, by in cls.PLANTED: + changed[row, 0] += np.float32(by) + tiles["ne.tif"] = piece(source, 0, 11, 10, 21, array=changed) + write_tiles(tmp_path / "disagree", tiles) + return tmp_path / "disagree" + + @staticmethod + def utm33(local: list[tuple[float, float]]) -> Ring: + return [(X0 + u, Y0 + v) for u, v in local] + + @staticmethod + def needed_nodes(polygon: Polygon, grid: Any) -> np.ndarray: + """The oracle's mask: `grid`'s nodes the domain grown by `grid`'s cell + diagonal, mitred, covers (closed).""" + grown = polygon.buffer(math.hypot(grid.delta_x, grid.delta_y), join_style="mitre") + r, c = np.indices((grid.rows, grid.cols)) + xs, ys = grid.x_min + c * grid.delta_x, grid.y_max - r * grid.delta_y + return np.asarray(shapely.covers(grown, shapely.points(xs, ys))) + + def opened(self, di: ModuleType, dm: ModuleType, tmp_path: Path, local: Any) -> Any: + dem = self.disagreeing(tmp_path) + return di.open_dem( + request(di, dem, domain=read(dm, tmp_path, self.utm33(local), "EPSG:25833")) + ) + + def bbox_of(self, di: ModuleType, mz: ModuleType, opened: Any, tmp_path: Path) -> Any: + x_min, y_min, x_max, y_max = opened.domain.polygon.bounds + bounds = mz.Bounds(x_min=x_min, y_min=y_min, x_max=x_max, y_max=y_max) + return di.open_dem(request(di, tmp_path / "disagree", bounds=bounds)) + + def tiles(self, tmp_path: Path) -> dict[str, Any]: + return {p.name: decoded(p) for p in sorted((tmp_path / "disagree").iterdir())} + + def test_a_seam_wholly_outside_the_needed_region_is_not_reported( + self, di: ModuleType, dm: ModuleType, tmp_path: Path + ) -> None: + """Ruling point 1: the planted overlap is inside the plan's rectangle, + both tiles are in the plan, and the region misses every planted node.""" + opened = self.opened(di, dm, tmp_path, self.STRIP) + m = opened.plan.meta + assert {"ne.tif", "nw.tif"} <= {t.name for t in opened.plan.tiles} + assert (m.x_min, m.y_max, m.rows, m.cols) == (X0, Y0, 21, 21) # the planted nodes in it + mask = self.needed_nodes(opened.domain.polygon, m) + assert not mask[1:6, 10].any() # the region misses them + assert mask[10, 10] # and crosses the column where the tiles agree + assert opened.seams == () + assert seams_of(self.tiles(tmp_path), m, mask) == [] + + def test_a_seam_the_needed_region_crosses_counts_only_its_nodes( + self, di: ModuleType, dm: ModuleType, tmp_path: Path + ) -> None: + """Ruling point 2: the L's region takes (100, -10) and (100, -20) and + leaves the other planted nodes, which the plan's rectangle holds.""" + opened = self.opened(di, dm, tmp_path, self.ELL) + m = opened.plan.meta + assert (m.x_min, m.y_max, m.rows, m.cols) == (X0, Y0, 21, 21) + mask = self.needed_nodes(opened.domain.polygon, m) + assert mask[1:6, 10].tolist() == [True, True, False, False, False] + assert seams_of(self.tiles(tmp_path), m, mask) == self.NORTH_TWO # the oracle + assert seams(opened) == self.NORTH_TWO + + def test_a_node_outside_the_polygon_inside_one_diagonal_counts( + self, di: ModuleType, dm: ModuleType, tmp_path: Path + ) -> None: + """The region is the grown polygon, not the polygon: every counted node + here is 7 m east of the domain, and bilinear z reads it there.""" + opened = self.opened(di, dm, tmp_path, self.WEST) + m = opened.plan.meta + polygon = opened.domain.polygon + for y in (Y0 - 10, Y0 - 20): + node = shapely.Point(X0 + 100, y) + assert not polygon.covers(node) + assert 0 < polygon.distance(node) < math.hypot(m.delta_x, m.delta_y) + assert seams_of(self.tiles(tmp_path), m, self.needed_nodes(polygon, m)) == self.NORTH_TWO + assert seams(opened) == self.NORTH_TWO + + @pytest.mark.parametrize("shape", ["STRIP", "ELL"]) + def test_without_a_domain_the_same_tiles_report_every_node( + self, di: ModuleType, dm: ModuleType, mz: ModuleType, tmp_path: Path, shape: str + ) -> None: + """Ruling: without `--domain` the report is unchanged. `--bbox` at the + domain's bounds has the same plan and counts all four planted nodes.""" + opened = self.opened(di, dm, tmp_path, getattr(self, shape)) + boxed = self.bbox_of(di, mz, opened, tmp_path) + assert boxed.plan == opened.plan + assert seams(boxed) == self.ALL + assert seams_of(self.tiles(tmp_path), boxed.tile.meta) == self.ALL + + @pytest.mark.parametrize("shape", ["STRIP", "ELL", "WEST"]) + def test_which_value_a_node_takes_is_unchanged_by_the_domain( + self, di: ModuleType, dm: ModuleType, mz: ModuleType, tmp_path: Path, shape: str + ) -> None: + """Ruling: the mask is the report's only. On the same plan the domain's + mosaic is `--bbox`'s bit for bit, and the midline oracle's.""" + opened = self.opened(di, dm, tmp_path, getattr(self, shape)) + boxed = self.bbox_of(di, mz, opened, tmp_path) + assert boxed.plan == opened.plan + assert same_array(np.asarray(opened.tile.array), np.asarray(boxed.tile.array)) + oracle = deepest_interior(self.tiles(tmp_path), opened.tile.meta) + assert np.array_equal(np.asarray(opened.tile.array, dtype=np.float64), oracle) diff --git a/tests/python/test_domain.py b/tests/python/test_domain.py index 74791e6a..1400799e 100644 --- a/tests/python/test_domain.py +++ b/tests/python/test_domain.py @@ -1,16 +1,25 @@ -"""`tin_engine.domain`: reading the domain polygon. Increment 16, R1. +"""`tin_engine.domain`: reading the domain polygon, and moving it into the DEM's CRS. -`docs/increments/16-domain-polygon.md`, R1, U1 (a) and U4 (a). The design -names `DomainPolygon` and `read_domain(path, crs)`; the extent check needs the -DEM's node rectangle, so the signature chosen here is:: +Increment 16, R1, U1 (a) and U4 (a), as changed by increment 15b +(`docs/increments/15-dem-mosaic.md` R9): the domain keeps its own CRS, and +`DomainPolygon.to_crs` replaces 16's must-match rule. Signatures, as pinned by +the 15b red suite ("Pinned by the red suite (15b)"):: - read_domain(path: Path, meta: RasterMeta, crs: str | None = None) -> DomainPolygon + read_domain(path: Path, crs: str | None = None) -> DomainPolygon + DomainPolygon.to_crs(dst: str | pyproj.CRS) -> DomainPolygon + check_extent(domain: DomainPolygon, meta: RasterMeta) -> None + +`read_domain` no longer takes the DEM: the domain is read before the tiles are +planned, because its bounds in the DEM's CRS choose them (R6). The extent +check runs after the transform, in `check_extent`, against the mosaic's node +rectangle (and the coverage check, in `open_dem`). `DomainPolygon` is frozen with `.polygon` (a shapely `Polygon`, oriented outer -counter-clockwise and holes clockwise) and `.epsg` (int). Every refusal raises -`DomainError`, a `ValueError` subclass (name chosen here), whose message the -CLI shows. `crs` is the `--domain-crs` text, `EPSG:n`; it is what a `.wkt` -file needs. +counter-clockwise and holes clockwise) and `.crs` (`str`, text pyproj parses to +the domain's CRS; it was `.epsg: int` in 16). Every refusal raises +`DomainError`, a `ValueError` subclass, whose message the CLI shows. `crs` is +the `--domain-crs` text, anything pyproj parses; it is what a `.wkt` file +needs, and for GeoJSON it must agree with the file's own (16, unchanged). The module is fetched inside a fixture, so its absence fails these tests and leaves the rest of the session collecting. @@ -24,13 +33,19 @@ from types import ModuleType from typing import Any +import numpy as np import pytest import shapely +from pyproj import CRS, Transformer from shapely.geometry import Polygon, mapping from tin_engine.io.models import RasterMeta UTM33 = "urn:ogc:def:crs:EPSG::25833" +# UTM 33 with the easting axis pointing west: `to_epsg()` still answers 25833 +# (at pyproj's default confidence), but it is not EPSG:25833, and the numbers +# differ in sign. A same-CRS test by EPSG code would skip its transform. +UTM33_WEST = "+proj=utm +zone=33 +ellps=GRS80 +units=m +axis=wnu +no_defs" # Nodes at x 500 000 .. 500 080 and y 6 599 970 .. 6 600 000. META = RasterMeta( @@ -95,19 +110,45 @@ def square(holes: bool = True) -> dict[str, Any]: def refused(domain: ModuleType, path: Path, *says: str, crs: str | None = None) -> None: with pytest.raises(domain.DomainError) as info: - domain.read_domain(path, META, crs) + domain.read_domain(path, crs) assert isinstance(info.value, ValueError) for word in says: assert word in str(info.value), str(info.value) +def transformed(src: str, dst: str, ring: list[tuple[float, float]]) -> list[tuple[float, float]]: + """pyproj's own `always_xy` transform of `ring`, vertex by vertex: the oracle.""" + t = Transformer.from_crs(src, dst, always_xy=True) + xs, ys = t.transform(np.array([p[0] for p in ring]), np.array([p[1] for p in ring])) + return [(float(x), float(y)) for x, y in zip(xs, ys, strict=True)] + + +def ring_set(ring: Any) -> set[tuple[float, float]]: + return {(float(x), float(y)) for x, y in ring.coords} + + +def the_crs(out: Any) -> CRS: + assert isinstance(out.crs, str) + return CRS.from_user_input(out.crs) + + +@pytest.fixture +def no_transformer(monkeypatch: pytest.MonkeyPatch) -> None: + """Any `Transformer.from_crs` from here on fails the test.""" + + def refuse(*args: Any, **kwargs: Any) -> Any: + raise AssertionError(f"Transformer.from_crs{args} was called") + + monkeypatch.setattr(Transformer, "from_crs", staticmethod(refuse)) + + class TestReading: @pytest.mark.parametrize("wrap", ["geometry", "feature", "collection"]) def test_geojson_as_geometry_feature_or_one_feature_collection( self, domain: ModuleType, tmp_path: Path, wrap: str ) -> None: - out = domain.read_domain(write(tmp_path, geojson(square(), wrap=wrap)), META) - assert out.epsg == 25833 + out = domain.read_domain(write(tmp_path, geojson(square(), wrap=wrap))) + assert the_crs(out) == CRS.from_epsg(25833) assert isinstance(out.polygon, Polygon) assert len(out.polygon.interiors) == 1 @@ -115,7 +156,8 @@ def test_geojson_as_geometry_feature_or_one_feature_collection( def test_both_spellings_of_the_crs_member( self, domain: ModuleType, tmp_path: Path, crs: str ) -> None: - assert domain.read_domain(write(tmp_path, geojson(square(), crs=crs)), META).epsg == 25833 + out = domain.read_domain(write(tmp_path, geojson(square(), crs=crs))) + assert the_crs(out) == CRS.from_epsg(25833) @pytest.mark.parametrize("flag", ["EPSG:25833", "epsg:25833"]) @pytest.mark.parametrize("member", [UTM33, "EPSG:25833"]) @@ -123,21 +165,55 @@ def test_geojson_with_an_agreeing_domain_crs( self, domain: ModuleType, tmp_path: Path, member: str, flag: str ) -> None: path = write(tmp_path, geojson(square(), crs=member)) - assert domain.read_domain(path, META, flag).epsg == 25833 + assert the_crs(domain.read_domain(path, flag)) == CRS.from_epsg(25833) def test_the_json_suffix_is_geojson(self, domain: ModuleType, tmp_path: Path) -> None: path = write(tmp_path, geojson(square()), name="d.json") - assert domain.read_domain(path, META).epsg == 25833 + assert the_crs(domain.read_domain(path)) == CRS.from_epsg(25833) def test_wkt_with_a_crs(self, domain: ModuleType, tmp_path: Path) -> None: path = write(tmp_path, Polygon(OUTER, [HOLE]).wkt, name="d.wkt") - out = domain.read_domain(path, META, "EPSG:25833") - assert out.epsg == 25833 + out = domain.read_domain(path, "EPSG:25833") + assert the_crs(out) == CRS.from_epsg(25833) assert len(out.polygon.interiors) == 1 + def test_geojson_without_crs_is_wgs84_by_the_standard( + self, domain: ModuleType, tmp_path: Path + ) -> None: + """RFC 7946: no `crs` member is WGS 84. 15b reads it as such, where 16 + refused it as a mismatch.""" + lon_lat = transformed("EPSG:25833", "EPSG:4326", OUTER) + out = domain.read_domain(write(tmp_path, geojson(dict(mapping(Polygon(lon_lat))), None))) + assert the_crs(out) == CRS.from_epsg(4326) + assert ring_set(out.polygon.exterior) == set(lon_lat) + + def test_another_utm_zone_is_read_in_its_own_crs( + self, domain: ModuleType, tmp_path: Path + ) -> None: + out = domain.read_domain(write(tmp_path, geojson(square(), crs="EPSG:25832"))) + assert the_crs(out) == CRS.from_epsg(25832) + + def test_wkt_in_a_geographic_crs(self, domain: ModuleType, tmp_path: Path) -> None: + lon_lat = transformed("EPSG:25833", "EPSG:4326", OUTER) + out = domain.read_domain(write(tmp_path, Polygon(lon_lat).wkt, name="d.wkt"), "EPSG:4326") + assert the_crs(out) == CRS.from_epsg(4326) + + def test_a_crs_without_an_epsg_code(self, domain: ModuleType, tmp_path: Path) -> None: + """R9: `crs: str`, "since a domain may have no EPSG code". 16 refused it.""" + lcc = "+proj=lcc +lat_1=60 +lat_2=65 +lat_0=62 +lon_0=15 +ellps=GRS80 +units=m +no_defs" + out = domain.read_domain(write(tmp_path, Polygon(OUTER).wkt, name="d.wkt"), lcc) + assert the_crs(out) == CRS.from_user_input(lcc) + + def test_ogc_crs84_member(self, domain: ModuleType, tmp_path: Path) -> None: + """Older GeoJSON writers name longitude-latitude WGS 84 this way.""" + lon_lat = transformed("EPSG:25833", "EPSG:4326", OUTER) + doc = geojson(dict(mapping(Polygon(lon_lat))), "urn:ogc:def:crs:OGC:1.3:CRS84") + out = domain.read_domain(write(tmp_path, doc)) + assert the_crs(out) == CRS.from_user_input("OGC:CRS84") + def test_vertices_are_used_as_given(self, domain: ModuleType, tmp_path: Path) -> None: """Directions 3 and 6: no densifying, simplifying or snapping.""" - out = domain.read_domain(write(tmp_path, geojson(square())), META) + out = domain.read_domain(write(tmp_path, geojson(square()))) assert set(out.polygon.exterior.coords) == set(OUTER) assert set(out.polygon.interiors[0].coords) == set(HOLE) assert len(out.polygon.exterior.coords) == len(OUTER) + 1 @@ -148,9 +224,151 @@ def test_orientation_is_outer_ccw_and_holes_cw( ) -> None: outer, hole = (OUTER[::-1], HOLE[::-1]) if reverse else (OUTER, HOLE) doc = geojson(dict(mapping(Polygon(outer, [hole])))) - out = domain.read_domain(write(tmp_path, doc), META) + out = domain.read_domain(write(tmp_path, doc)) + assert out.polygon.exterior.is_ccw + assert not out.polygon.interiors[0].is_ccw + + def test_a_vertex_outside_any_dem_is_not_reading_s_business( + self, domain: ModuleType, tmp_path: Path + ) -> None: + """The extent is checked after the transform, against the tiles (R9).""" + far = [(x + 1e6, y) for x, y in OUTER] + out = domain.read_domain(write(tmp_path, geojson(dict(mapping(Polygon(far)))))) + assert ring_set(out.polygon.exterior) == set(far) + + def test_the_model_is_frozen(self, domain: ModuleType, tmp_path: Path) -> None: + out = domain.read_domain(write(tmp_path, geojson(square()))) + # Pydantic's ValidationError is a ValueError; a frozen dataclass's + # FrozenInstanceError is an AttributeError. + with pytest.raises((ValueError, AttributeError, TypeError)): + out.crs = "EPSG:4326" + + def test_the_must_match_rule_is_gone(self, domain: ModuleType) -> None: + """R9: `check_crs` is replaced by the transform, and `epsg: int` by `crs: str`.""" + assert not hasattr(domain, "check_crs") + assert "epsg" not in domain.DomainPolygon.model_fields + assert "crs" in domain.DomainPolygon.model_fields + + +class TestToCrs: + """R9: vertices only, each through pyproj's `always_xy` transform, exactly once.""" + + def lon_lat_domain(self, domain: ModuleType, tmp_path: Path) -> Any: + outer = transformed("EPSG:25833", "EPSG:4326", OUTER) + hole = transformed("EPSG:25833", "EPSG:4326", HOLE) + doc = geojson(dict(mapping(Polygon(outer, [hole]))), crs=None) + return domain.read_domain(write(tmp_path, doc)) + + def test_wgs84_to_utm33_is_pyprojs_transform_of_every_vertex( + self, domain: ModuleType, tmp_path: Path + ) -> None: + read = self.lon_lat_domain(domain, tmp_path) + out = read.to_crs("EPSG:25833") + assert the_crs(out) == CRS.from_epsg(25833) + for before, after in ( + (read.polygon.exterior, out.polygon.exterior), + (read.polygon.interiors[0], out.polygon.interiors[0]), + ): + expected = transformed("EPSG:4326", "EPSG:25833", list(before.coords)) + assert ring_set(after) == set(expected) + # Vertices only: not densified, nothing dropped. + assert len(after.coords) == len(before.coords) + + def test_the_result_lands_where_the_utm_polygon_was( + self, domain: ModuleType, tmp_path: Path + ) -> None: + """A round trip through longitude and latitude is not exact, but it is + within a micrometre; a swapped axis would be thousands of km off.""" + out = self.lon_lat_domain(domain, tmp_path).to_crs("EPSG:25833") + got = sorted(ring_set(out.polygon.exterior)) + for (x, y), (ex, ey) in zip(got, sorted(set(OUTER)), strict=True): + assert x == pytest.approx(ex, abs=1e-6) + assert y == pytest.approx(ey, abs=1e-6) + + def test_the_source_is_unchanged(self, domain: ModuleType, tmp_path: Path) -> None: + read = self.lon_lat_domain(domain, tmp_path) + before = list(read.polygon.exterior.coords) + read.to_crs("EPSG:25833") + assert list(read.polygon.exterior.coords) == before + assert the_crs(read) == CRS.from_epsg(4326) + + @pytest.mark.parametrize("src", ["EPSG:25832", "EPSG:3035"]) + def test_projected_to_projected(self, domain: ModuleType, tmp_path: Path, src: str) -> None: + ring = transformed("EPSG:25833", src, OUTER) + read = domain.read_domain(write(tmp_path, Polygon(ring).wkt, name="d.wkt"), src) + out = read.to_crs("EPSG:25833") + assert ring_set(out.polygon.exterior) == set(transformed(src, "EPSG:25833", ring)) + + def test_a_crs_object_as_the_destination(self, domain: ModuleType, tmp_path: Path) -> None: + read = self.lon_lat_domain(domain, tmp_path) + by_text = read.to_crs("EPSG:25833") + by_object = read.to_crs(CRS.from_epsg(25833)) + assert list(by_object.polygon.exterior.coords) == list(by_text.polygon.exterior.coords) + + def test_the_winding_contract_survives_a_mirroring_transform( + self, domain: ModuleType, tmp_path: Path + ) -> None: + """From a west-pointing easting, the transform mirrors the ring; the + result is re-oriented, outer counter-clockwise and holes clockwise.""" + outer = [(-x, y) for x, y in OUTER] + hole = [(-x, y) for x, y in HOLE] + path = write(tmp_path, Polygon(outer, [hole]).wkt, name="d.wkt") + out = domain.read_domain(path, UTM33_WEST).to_crs("EPSG:25833") assert out.polygon.exterior.is_ccw assert not out.polygon.interiors[0].is_ccw + assert ring_set(out.polygon.exterior) == set(transformed(UTM33_WEST, "EPSG:25833", outer)) + + def test_the_same_epsg_code_is_not_the_same_crs( + self, domain: ModuleType, tmp_path: Path + ) -> None: + """Same-CRS is decided by CRS equality, not by `to_epsg()`: this CRS's + `to_epsg()` is 25833, and skipping its transform would put the domain + a thousand kilometres west of where it is.""" + assert CRS.from_user_input(UTM33_WEST).to_epsg() == 25833 # the trap is real + outer = [(-x, y) for x, y in OUTER] + out = domain.read_domain(write(tmp_path, Polygon(outer).wkt, name="d.wkt"), UTM33_WEST) + moved = out.to_crs("EPSG:25833") + assert ring_set(moved.polygon.exterior) == set(transformed(UTM33_WEST, "EPSG:25833", outer)) + assert all(x > 0 for x, _ in moved.polygon.exterior.coords) + + @pytest.mark.parametrize("dst", ["EPSG:25833", "urn:ogc:def:crs:EPSG::25833"]) + def test_the_same_crs_is_not_transformed( + self, domain: ModuleType, tmp_path: Path, no_transformer: None, dst: str + ) -> None: + """A domain already in the DEM's CRS keeps its coordinates bit for bit, + and no transformer is made (the mesh must be 16's, bit for bit).""" + read = domain.read_domain(write(tmp_path, geojson(square(), crs=UTM33))) + out = read.to_crs(dst) + assert list(out.polygon.exterior.coords) == list(read.polygon.exterior.coords) + assert list(out.polygon.interiors[0].coords) == list(read.polygon.interiors[0].coords) + assert the_crs(out) == CRS.from_epsg(25833) + + def test_a_vertex_with_no_image_is_refused_naming_both_crss( + self, domain: ModuleType, tmp_path: Path + ) -> None: + """Latitude 95 has no image in UTM 33 (pyproj answers inf). Refused, + at reading or at the transform, never passed on as a coordinate.""" + ring = [(14.9, 59.4), (15.1, 59.4), (15.1, 95.0), (14.9, 59.6)] + path = write(tmp_path, geojson(dict(mapping(Polygon(ring))), crs=None)) + with pytest.raises(domain.DomainError) as info: + domain.read_domain(path).to_crs("EPSG:25833") + for word in ("4326", "25833"): + assert word in str(info.value), info.value + + def test_utm_numbers_in_a_file_without_a_crs_member_are_refused( + self, domain: ModuleType, tmp_path: Path + ) -> None: + """The commonest wrong file: UTM coordinates, no `crs` member, so WGS 84 + by the standard. Its "longitudes" are 500 000.""" + path = write(tmp_path, geojson(square(), crs=None)) + with pytest.raises(domain.DomainError) as info: + domain.read_domain(path).to_crs("EPSG:25833") + for word in ("4326", "25833"): + assert word in str(info.value), info.value + + +class TestCheckExtent: + """U4 (a), now after the transform, against the mosaic's node rectangle.""" def test_a_vertex_on_the_border_of_the_node_rectangle_is_inside( self, domain: ModuleType, tmp_path: Path @@ -161,15 +379,30 @@ def test_a_vertex_on_the_border_of_the_node_rectangle_is_inside( (500_080.0, 6_600_000.0), (500_000.0, 6_600_000.0), ] - out = domain.read_domain(write(tmp_path, geojson(dict(mapping(Polygon(corners))))), META) - assert out.polygon.bounds == (500_000.0, 6_599_970.0, 500_080.0, 6_600_000.0) + out = domain.read_domain(write(tmp_path, geojson(dict(mapping(Polygon(corners)))))) + assert domain.check_extent(out, META) is None - def test_the_model_is_frozen(self, domain: ModuleType, tmp_path: Path) -> None: - out = domain.read_domain(write(tmp_path, geojson(square())), META) - # Pydantic's ValidationError is a ValueError; a frozen dataclass's - # FrozenInstanceError is an AttributeError. - with pytest.raises((ValueError, AttributeError, TypeError)): - out.epsg = 4326 + @pytest.mark.parametrize("dx", [1.0, 1e-6]) + def test_a_vertex_outside_the_node_rectangle( + self, domain: ModuleType, tmp_path: Path, dx: float + ) -> None: + """Refused, naming the vertex and saying it is outside.""" + outer = [*OUTER[:1], (500_080.0 + dx, 6_599_975.0), *OUTER[2:]] + out = domain.read_domain(write(tmp_path, geojson(dict(mapping(Polygon(outer)))))) + with pytest.raises(domain.DomainError) as info: + domain.check_extent(out, META) + assert "outside" in str(info.value) + assert "500080" in str(info.value).replace(",", "").replace(" ", "") + + def test_a_hole_vertex_is_checked_too(self, domain: ModuleType, tmp_path: Path) -> None: + corners = [ + (500_000.0, 6_599_970.0), + (500_080.0, 6_599_970.0), + (500_080.0, 6_600_000.0), + (500_000.0, 6_600_000.0), + ] + out = domain.read_domain(write(tmp_path, geojson(dict(mapping(Polygon(corners, [HOLE])))))) + assert domain.check_extent(out, META) is None class TestRefusals: @@ -205,14 +438,6 @@ def test_invalid_names_the_reason(self, domain: ModuleType, tmp_path: Path) -> N doc = geojson({"type": "Polygon", "coordinates": [[*bowtie, bowtie[0]]]}) refused(domain, write(tmp_path, doc), "Self-intersection") - def test_geojson_without_crs_is_wgs84_by_the_standard( - self, domain: ModuleType, tmp_path: Path - ) -> None: - refused(domain, write(tmp_path, geojson(square(), crs=None)), "4326", "25833") - - def test_a_mismatched_epsg(self, domain: ModuleType, tmp_path: Path) -> None: - refused(domain, write(tmp_path, geojson(square(), crs="EPSG:25832")), "25832", "25833") - @pytest.mark.parametrize( ("member", "flag", "says"), [ @@ -229,7 +454,8 @@ def test_geojson_with_a_disagreeing_domain_crs( flag: str, says: tuple[str, ...], ) -> None: - """The flag never overrides the file's own CRS: a disagreement is refused. + """The flag never overrides the file's own CRS: a disagreement is refused + (16, unchanged by 15b: it is the file against the flag, not the DEM). The third row is a file with no ``crs`` member, which RFC 7946 makes WGS 84; ``--domain-crs`` does not supply one.""" @@ -238,16 +464,14 @@ def test_geojson_with_a_disagreeing_domain_crs( def test_wkt_without_a_crs(self, domain: ModuleType, tmp_path: Path) -> None: refused(domain, write(tmp_path, Polygon(OUTER).wkt, name="d.wkt")) - def test_wkt_with_a_mismatched_crs(self, domain: ModuleType, tmp_path: Path) -> None: - refused(domain, write(tmp_path, Polygon(OUTER).wkt, name="d.wkt"), "4326", crs="EPSG:4326") - - @pytest.mark.parametrize("dx", [1.0, 1e-6]) - def test_a_vertex_outside_the_node_rectangle( - self, domain: ModuleType, tmp_path: Path, dx: float + @pytest.mark.parametrize("bad", ["EPSG:999999", "not a crs"]) + def test_an_unknown_domain_crs_is_named( + self, domain: ModuleType, tmp_path: Path, bad: str ) -> None: - """U4 (a): refused, naming the vertex and saying it is outside.""" - outer = [*OUTER[:1], (500_080.0 + dx, 6_599_975.0), *OUTER[2:]] - refused(domain, write(tmp_path, geojson(dict(mapping(Polygon(outer))))), "outside") + refused(domain, write(tmp_path, Polygon(OUTER).wkt, name="d.wkt"), bad, crs=bad) + + def test_an_unknown_crs_member_is_named(self, domain: ModuleType, tmp_path: Path) -> None: + refused(domain, write(tmp_path, geojson(square(), crs="EPSG:999999")), "EPSG:999999") def test_an_unknown_suffix(self, domain: ModuleType, tmp_path: Path) -> None: refused(domain, write(tmp_path, geojson(square()), name="d.shp"), ".shp") diff --git a/tests/python/test_refine_golden.py b/tests/python/test_refine_golden.py index e9ff2919..6136d4a5 100644 --- a/tests/python/test_refine_golden.py +++ b/tests/python/test_refine_golden.py @@ -74,7 +74,7 @@ def refined( chains = [(ring, ChainRole.Outer, 0)] else: path = geojson(tmp_path / "quarter.geojson", quarter_circle()) - xy, chains, _ = _domain_chains(read_domain(path, tile.meta), path.name) + xy, chains, _ = _domain_chains(read_domain(path), path.name) run = _engine(xy, chains, True, DEFAULT_SNAP_SPACING) assert run.mesh is not None and run.noded is not None, run.message edges, masks = _constraint_arrays(run.mesh, run.noded)