From d97fcc7ad42eaf21a456320d6b0eef16751bafaf Mon Sep 17 00:00:00 2001 From: Zwitschi Date: Fri, 6 Nov 2026 23:01:21 +0100 Subject: [PATCH] geo: acquire the 3DEP 1 m DEM window in EPSG:26911 550x350 m clipped to the block, already in the project CRS, 984 KB, committed through LFS. Verified against the service's own point query and the earlier probe: 115.911/119.422/117.593 m, agreement within 6 cm. Adds Geo/scripts/dem_stats.py, a stdlib tiled-TIFF reader, because GDAL cannot run on the host and the exported raster is tiled rather than strip-organised. --- Geo/README.md | 2 +- Geo/manifests/2026-11-06-3dep-dem.md | 43 ++++++ Geo/raw/dem-3dep-1m-26911-2026-11-06.tif | 3 + Geo/scripts/dem_stats.py | 173 +++++++++++++++++++++++ Geo/verification/2026-11-06-elevation.md | 13 +- ROADMAP.md | 4 +- THIRD_PARTY_NOTICES.md | 2 +- 7 files changed, 233 insertions(+), 7 deletions(-) create mode 100644 Geo/manifests/2026-11-06-3dep-dem.md create mode 100644 Geo/raw/dem-3dep-1m-26911-2026-11-06.tif create mode 100755 Geo/scripts/dem_stats.py diff --git a/Geo/README.md b/Geo/README.md index 0bb9b64..c4c396d 100644 --- a/Geo/README.md +++ b/Geo/README.md @@ -13,7 +13,7 @@ OSM is layout data only. It is not authoritative for façade detail, height, or | `Geo/manifests/` | One manifest per saved extract (see below) | | `Geo/scripts/` | Repeatable fetch and conversion helpers | -Current state: the extract, both El Capitan footprints, the demo bounds, the building sets, the street profile, and the designated crossing are recorded in [`verification/`](verification/). The local CRS and origin are fixed in [`manifests/local-origin.md`](manifests/local-origin.md), with the projected easting/northing still to be computed. Elevation is measured (3.94 m relief, 3DEP chosen) but the DEM tile is not downloaded; the hero façade width remains open. +Current state: the extract, both El Capitan footprints, the demo bounds, the building sets, the street profile, and the designated crossing are recorded in [`verification/`](verification/). The local CRS and origin are fixed in [`manifests/local-origin.md`](manifests/local-origin.md), with the projected easting/northing still to be computed. Elevation is measured and the 3DEP 1 m DEM window is committed (EPSG:26911, NAVD88, 984 KB). The hero façade width remains open. ## Required manifest diff --git a/Geo/manifests/2026-11-06-3dep-dem.md b/Geo/manifests/2026-11-06-3dep-dem.md new file mode 100644 index 0000000..659cabf --- /dev/null +++ b/Geo/manifests/2026-11-06-3dep-dem.md @@ -0,0 +1,43 @@ +# Elevation source — USGS 3DEP 1 m DEM, block window + +| Field | Value | +| --- | --- | +| Source | USGS 3DEP Elevation ImageServer, `elevation.nationalmap.gov`, `/arcgis/rest/services/3DEPElevation/ImageServer/exportImage` | +| Product | 3DEP seamless 1 m (CONUS), raster id 129943 | +| Fetched | 2026-11-06 22:59 CET | +| Request | `f=image&bbox=376100,3774050,376650,3774400&bboxSR=26911&imageSR=26911&size=550,350&format=tiff&pixelType=F32&noData=-9999&interpolation=RSP_BilinearInterpolation` | +| Bounds (EPSG:26911) | E 376 100 – 376 650, N 3 774 050 – 3 774 400 | +| Bounds (WGS84, approx.) | lat 34.1003 – 34.1032, lon −118.3426 – −118.3370 | +| Grid | 550 × 350 px, 1 m pixel, float32, tiled TIFF, nodata −9999 | +| File | `Geo/raw/dem-3dep-1m-26911-2026-11-06.tif`, 984 338 bytes | +| SHA-256 | `e571c14c7230d9048697db2893095e1700468e45d4fccf258ceb4d4e9291bdc5` | +| Horizontal CRS | EPSG:26911 requested and returned (GeoKey 3072 = 26911, model type projected) | +| Vertical datum | NAVD88 orthometric — 3DEP CONUS convention; the service does not restate it in the file, so this is a documented assumption | +| Licence | US Government work, public domain | + +The window is the block plus a ~100 m margin on the east and ~70 m on the other three sides, which is enough to carry the grade through the junctions at both ends without dragging in the hills north of the boulevard. + +## Sampling cross-check + +The raster is reprojected by the service, not by us, so it was checked against the service's own point query (`getSamples`) and against the independent point service used in `../verification/2026-11-06-elevation.md`: + +| Point | Raster pixel | Service `getSamples` | Earlier point probe | +| --- | ---: | ---: | ---: | +| Hollywood × Highland (origin) | 115.911 | 115.910 | 115.911 | +| Hollywood × Orange Dr | 119.422 | 119.418 | 119.418 | +| Mid-block crossing | 117.593 | — | 117.53 | + +Agreement is within 4 cm on the raster side and 6 cm against the point probe, so the file's georeferencing and signage are consistent. + +## What the window contains + +Min 110.46 m, max 127.73 m, mean 117.43 m over the 550 × 350 m window — a 17.28 m spread. That spread is **not** the street relief: it includes building tops and the rising ground north of the boulevard. The walkable relief measured along the street is 3.94 m (`../verification/2026-11-06-elevation.md`), and that is the number the demo's ground mesh must honour. + +## Reading the file + +`Geo/scripts/dem_stats.py` parses it with the standard library. GDAL cannot run on the host (its flatpak libraries do not resolve outside the sandbox), and the exported TIFF is **tiled**, not strip-organised — a strip reader silently scrambles the rows. If another tool misreads it, suspect the tile layout before the data. + +## Open + +- Reprojecting is not needed: the service already returned EPSG:26911. Anything that samples this file must assume NAVD88 and treat −9999 as nodata. +- The window stops 30 m north of the Chinese Theatre's roof line; if the background massing needs ground beyond it, extend the bbox rather than extrapolating. diff --git a/Geo/raw/dem-3dep-1m-26911-2026-11-06.tif b/Geo/raw/dem-3dep-1m-26911-2026-11-06.tif new file mode 100644 index 0000000..d2b26c4 --- /dev/null +++ b/Geo/raw/dem-3dep-1m-26911-2026-11-06.tif @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:e571c14c7230d9048697db2893095e1700468e45d4fccf258ceb4d4e9291bdc5 +size 984338 diff --git a/Geo/scripts/dem_stats.py b/Geo/scripts/dem_stats.py new file mode 100755 index 0000000..6796e92 --- /dev/null +++ b/Geo/scripts/dem_stats.py @@ -0,0 +1,173 @@ +#!/usr/bin/env python3 +"""Read a plain float32 GeoTIFF, report its georeferencing and elevation stats. + +Usage: + Geo/scripts/dem_stats.py Geo/raw/dem-3dep-1m-26911-2026-11-06.tif + Geo/scripts/dem_stats.py --sample 376515.389,3774224.083 + +Why this exists: GDAL is not runnable on the host (its flatpak libraries do not +resolve outside the sandbox), so this reads the strip-based GeoTIFF that the +3DEP export service produces with the standard library only. It handles +uncompressed, single-sample, strip-organised float32 images with the usual +GeoTIFF tags; it is a small reader, not a general GeoTIFF implementation. + +Values are orthometric metres (NAVD88 for 3DEP CONUS). Nodata is excluded from +the statistics. +""" + +import argparse +import struct +import sys + +TYPE_SIZES = {1: 1, 2: 1, 3: 2, 4: 4, 5: 8, 7: 1, 11: 4, 12: 8, 16: 8} + + +def read_tiff(path): + with open(path, "rb") as handle: + blob = handle.read() + + if blob[:2] == b"II": + endian = "<" + elif blob[:2] == b"MM": + endian = ">" + else: + raise SystemExit("not a TIFF file") + if struct.unpack(endian + "H", blob[2:4])[0] != 42: + raise SystemExit("not a classic TIFF (BigTIFF is unsupported)") + + ifd_offset = struct.unpack(endian + "I", blob[4:8])[0] + count = struct.unpack(endian + "H", blob[ifd_offset : ifd_offset + 2])[0] + entries = {} + for i in range(count): + base = ifd_offset + 2 + i * 12 + tag, typ, n = struct.unpack(endian + "HHI", blob[base : base + 8]) + size = TYPE_SIZES.get(typ, 1) * n + raw = blob[base + 8 : base + 8 + 4] if size <= 4 else blob[ + struct.unpack(endian + "I", blob[base + 8 : base + 12])[0] : + ][:size] + fmt = {1: "B", 3: "H", 4: "I", 7: "B", 11: "f", 12: "d", 2: "c"}.get(typ, "B") + if typ == 2: + values = raw[:size].split(b"\0")[0].decode("latin-1") + elif fmt == "c": + values = raw + else: + values = list(struct.unpack(endian + fmt * n, raw[: struct.calcsize(endian + fmt * n)])) + entries[tag] = values + + width = entries[256][0] + height = entries[257][0] + bits = entries.get(258, [8])[0] + compression = entries.get(259, [1])[0] + sample_format = entries.get(339, [1])[0] + if compression != 1 or bits != 32 or sample_format != 3: + raise SystemExit( + f"unsupported TIFF layout: compression={compression} bits={bits} sample_format={sample_format}" + ) + + if 322 in entries: + # tiled layout: tags 322/323 give tile size, 324/325 the offsets and byte counts + tile_w = entries[322][0] + tile_h = entries[323][0] + offsets = entries[324] + counts = entries[325] + across = (width + tile_w - 1) // tile_w + down = (height + tile_h - 1) // tile_h + if len(offsets) < across * down: + raise SystemExit(f"tiled TIFF with {len(offsets)} tiles, expected {across * down}") + grid = [[0.0] * width for _ in range(height)] + for index, (offset, nbytes) in enumerate(zip(offsets, counts)): + tile_row, tile_col = divmod(index, across) + tile = struct.unpack(endian + "f" * (tile_w * tile_h), blob[offset : offset + nbytes]) + for r in range(tile_h): + target_row = tile_row * tile_h + r + if target_row >= height: + break + base = r * tile_w + for c in range(tile_w): + target_col = tile_col * tile_w + c + if target_col >= width: + break + grid[target_row][target_col] = tile[base + c] + flat = [value for row in grid for value in row] + return entries, width, height, flat + + # strip layout + offsets = entries[273] + counts = entries[279] + rows_per_strip = entries.get(278, [height])[0] + chunks = [blob[offset : offset + nbytes] for offset, nbytes in zip(offsets, counts)] + pixels = b"".join(chunks) + values = list(struct.unpack(endian + "f" * (width * height), pixels[: width * height * 4])) + return entries, width, height, values + + +def geokeys(entries): + keys = {} + raw = entries.get(34735) + if not raw: + return keys + # GeoKeyDirectory: version, revision, minor, count, then 4 shorts per key + for i in range(4, len(raw), 4): + key, tag_location, count, value = raw[i : i + 4] + keys[key] = (tag_location, count, value) + return keys + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("tif") + ap.add_argument("--sample", action="append", default=[], help="E,N to sample, repeatable") + args = ap.parse_args() + + entries, width, height, values = read_tiff(args.tif) + scale = entries.get(33550) + tiepoint = entries.get(33922) + keys = geokeys(entries) + nodata = None + if 42113 in entries: + try: + nodata = float(entries[42113]) + except ValueError: + nodata = None + + print(f"size: {width} x {height} px, float32, uncompressed") + if nodata is not None: + print(f"nodata: {nodata}") + + valid = [v for v in values if nodata is None or v != nodata] + if valid: + print( + f"elevation: min {min(valid):.2f} m, max {max(valid):.2f} m, " + f"relief {max(valid) - min(valid):.2f} m, mean {sum(valid)/len(valid):.2f} m" + ) + print(f"valid pixels: {len(valid)} of {width*height}") + + if scale and tiepoint: + sx, sy = scale[0], scale[1] + px, py, pz, ex, ey = tiepoint[0], tiepoint[1], tiepoint[2], tiepoint[3], tiepoint[4] + x0 = ex - px * sx + y0 = ey + py * sy # north-up: tiepoint row 0 is the top edge + print(f"pixel size: {sx} x {sy} m") + print(f"extent: E {x0:.1f} .. {x0 + width*sx:.1f}, N {y0 - height*sy:.1f} .. {y0:.1f}") + crs = keys.get(3072, (None, None, None))[2] + model = keys.get(1024, (None, None, None))[2] + print(f"model type: {model} (1 = projected), projected CRS: EPSG:{crs}") + + for spec in args.sample: + try: + e, n = (float(x) for x in spec.split(",")) + except ValueError: + raise SystemExit(f"bad --sample value: {spec!r}") + col = int((e - x0) / sx) + row = int((y0 - n) / sy) + if not (0 <= col < width and 0 <= row < height): + print(f" sample E{e:.1f} N{n:.1f}: outside the raster") + continue + value = values[row * width + col] + print(f" sample E{e:.3f} N{n:.3f} -> row {row}, col {col}: {value:.3f} m") + else: + print("no pixel scale or tiepoint found; georeferencing missing", file=sys.stderr) + + +if __name__ == "__main__": + main() diff --git a/Geo/verification/2026-11-06-elevation.md b/Geo/verification/2026-11-06-elevation.md index 4dc2585..2e58390 100644 --- a/Geo/verification/2026-11-06-elevation.md +++ b/Geo/verification/2026-11-06-elevation.md @@ -42,8 +42,15 @@ Consequence for `0.1.0`: the demo cannot use a flat ground plane. The street and Points inside the hero footprint returned 114.53–117.13 m — a 2.6 m spread within one building outline, which is not ground reference. 3DEP can carry structures and voids there. Those four diagnostic points are excluded from the relief figure; use the street and sidewalk lines, not building interiors, for any ground surface. +## DEM acquired + +Downloaded 2026-11-06 from the 3DEP export service, already in EPSG:26911, clipped to E 376 100–376 650 / N 3 774 050–3 774 400 at 1 m. Manifest, checksum, and full request: [`../manifests/2026-11-06-3dep-dem.md`](../manifests/2026-11-06-3dep-dem.md). File is committed — it is 984 KB, public domain, and the blockout needs it. + +Cross-check at three points: raster 115.911 / 119.422 / 117.593 m against the service's own point query 115.910 / 119.418 and this document's probe 115.911 / 119.418 / 117.53. Agreement is centimetres, so the georeferencing holds. + +Window statistics: min 110.46 m, max 127.73 m, spread 17.28 m — building tops and the ground rising north of the boulevard. The **walkable relief stays 3.94 m**; the demo ground mesh follows the street, not the window envelope. + ## Next -1. Download the 3DEP tile(s) covering the demo bounds, keep the source file (size is larger than the OSM extract; decide commit vs local-only with a manifest and checksum at that point). -2. Reproject to EPSG:26911, clip to the block plus margin, and record the vertical datum and the clip bounds. -3. Sample the kerb line and the crossing ramps specifically — the grid above is a street-level average, not a kerb profile. +1. Sample the kerb line and the crossing ramps specifically — the grid above is a street-level average, not a kerb profile. +2. Build the ground mesh from the DEM along the street, not from a flat plane, and assert the 3.94 m fall survives the Blender→Unreal path. diff --git a/ROADMAP.md b/ROADMAP.md index 2d74036..c8f0dbc 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -27,8 +27,8 @@ Status detail in [Docs/setup.md](Docs/setup.md); workstation audit in [Docs/envi ### Geographic data - [x] Save a versioned OSM extract; document date, bounds, and source in `Geo/`. -- [-] Acquire separately licensed elevation data if needed. Needed: measured relief is 3.94 m over the block (`Geo/verification/2026-11-06-elevation.md`). Source chosen: USGS 3DEP 1 m DEM, public domain, vertical datum NAVD88. -- [ ] Download the 3DEP 1 m DEM for the demo bounds, reproject it to EPSG:26911 with the recorded origin and vertical datum, and record the source, clip bounds, and checksum in `Geo/`. +- [x] Acquire separately licensed elevation data if needed. Needed and acquired: USGS 3DEP 1 m DEM in EPSG:26911, public domain, NAVD88, 984 KB committed (`Geo/manifests/2026-11-06-3dep-dem.md`); measured relief 3.94 m. +- [x] Download the 3DEP 1 m DEM for the demo bounds, reproject it to EPSG:26911 with the recorded origin and vertical datum, and record the source, clip bounds, and checksum in `Geo/`. - [x] In QGIS, set a suitable local projected CRS (e.g. EPSG:26911) and record a local origin. EPSG:26911, origin and projected easting/northing are recorded in `Geo/manifests/local-origin.md` with PROJ 9.9.0 control points; QGIS 4.2.3 is installed for the DEM reprojection step. - [x] Verify the address, building footprints, and demo boundaries for 6840 Hollywood Blvd against source data. - [ ] Import OSM blockouts through Blosm; clean and detail them in Blender. diff --git a/THIRD_PARTY_NOTICES.md b/THIRD_PARTY_NOTICES.md index 6c930cd..04916cf 100644 --- a/THIRD_PARTY_NOTICES.md +++ b/THIRD_PARTY_NOTICES.md @@ -21,7 +21,7 @@ Status values: **verified** (licence confirmed and permits our use), **pending** | --- | --- | --- | --- | | OSM extract — Hollywood block, 2026-11-06 | OpenStreetMap API v0.6, `Geo/raw/osm-api-map-2026-11-06.osm` | ODbL 1.0 | verified — manifest and SHA-256 in `Geo/manifests/` | | OSM building ways 424024147, 424024148 | OpenStreetMap API v0.6, `Geo/raw/*.json` | ODbL 1.0 | verified | -| Elevation data — USGS 3DEP 1 m DEM | `epqs.nationalmap.gov` point service sampled; tile not yet downloaded | US Government work, public domain | verified as a source; acquisition and manifest pending | +| Elevation data — USGS 3DEP 1 m DEM | `elevation.nationalmap.gov` 3DEP ImageServer export, `Geo/raw/dem-3dep-1m-26911-2026-11-06.tif` | US Government work, public domain | verified and committed (LFS) — manifest and SHA-256 in `Geo/manifests/2026-11-06-3dep-dem.md` | Attribution string for any OSM-derived output: "© OpenStreetMap contributors", .