#!/usr/bin/env python3 """Prepare DTK25 Blätter into an Übersicht and optional Kacheln. Reads each Blatt's GeoTIFF tags with the standard library (no shell, no GDAL). Every Blatt must be 9600 by 9600 pixels at 1.25 m in EPSG:25832, spanning 12 km by 12 km, and must have a matching world file. The Übersicht lists every Blatt. Kacheln are written only when ``--tiles`` is passed, for one Blatt unless ``--all`` is given. ``--tiles public`` writes ``public/kacheln``. ``gdalwarp`` and ``gdal_translate`` are used when both are on PATH. Otherwise Pillow opens the LZW palette TIFF, keeps the native pixels at 1.25 m, and resizes each coarser step. Resolutions are metres per pixel, finest first. Each coarser step doubles the ground resolution, from the native sheet pixel up to 160 m. The phone's furthest ground span is app configuration and is not read here; changing that span does not require a new tile set. """ from __future__ import annotations import argparse import json import math import os import shutil import struct import subprocess import sys import tempfile from dataclasses import dataclass from pathlib import Path TILE_SIZE_PX = 256 NATIVE_RESOLUTION = 1.25 EXPECTED_PIXELS = 9600 EXPECTED_SPAN_M = 12_000 EXPECTED_EPSG = 25832 SCALE_TOLERANCE = 1e-6 SPAN_TOLERANCE_M = 0.01 PIXEL_TOLERANCE = 1e-4 RESOLUTIONS = [1.25, 2.5, 5, 10, 20, 40, 80, 160] TAG_WIDTH = 256 TAG_HEIGHT = 257 TAG_PIXEL_SCALE = 33550 TAG_TIEPOINT = 33922 TAG_GEOKEY = 34735 GEOKEY_PROJECTED_CS = 3072 TYPE_SIZE = { 1: 1, 2: 1, 3: 2, 4: 4, 5: 8, 6: 1, 7: 1, 8: 2, 9: 4, 10: 8, 11: 4, 12: 8, } @dataclass(frozen=True) class Blatt: path: Path min_e: float min_n: float max_e: float max_n: float def as_json(self) -> dict[str, int | float]: return { "minE": json_number(self.min_e), "minN": json_number(self.min_n), "maxE": json_number(self.max_e), "maxN": json_number(self.max_n), } def json_number(value: float) -> int | float: rounded = round(value) if abs(value - rounded) <= SCALE_TOLERANCE: return int(rounded) return round(value, 6) def snap(value: float, what: str) -> int: nearest = round(value) if abs(value - nearest) > PIXEL_TOLERANCE: raise SystemExit(f"{what}: {value} is not an integer pixel") return nearest def floor_index(value: float, size: float) -> int: quotient = value / size nearest = round(quotient) if abs(quotient - nearest) <= SCALE_TOLERANCE: return nearest return math.floor(quotient) def ceil_index(value: float, size: float) -> int: quotient = value / size nearest = round(quotient) if abs(quotient - nearest) <= SCALE_TOLERANCE: return nearest return math.ceil(quotient) def tile_indices(min_v: float, max_v: float, ground: float) -> range: """Columns or rows whose tiles intersect [min_v, max_v).""" start = floor_index(min_v, ground) stop = ceil_index(max_v, ground) - 1 return range(start, stop + 1) def intersects(min_e: float, min_n: float, max_e: float, max_n: float, blatt: Blatt) -> bool: return min_e < blatt.max_e and max_e > blatt.min_e and min_n < blatt.max_n and max_n > blatt.min_n def read_values(handle, endian: str, typ: int, count: int, raw: bytes) -> tuple: if typ not in TYPE_SIZE: raise ValueError(f"unsupported TIFF type {typ}") size = TYPE_SIZE[typ] * count if size <= 4: blob = raw[:size] else: (offset,) = struct.unpack(endian + "I", raw) handle.seek(offset) blob = handle.read(size) if len(blob) != size: raise ValueError(f"TIFF value at {offset} is short ({len(blob)} of {size} bytes)") if typ == 3: return struct.unpack(endian + f"{count}H", blob) if typ == 4: return struct.unpack(endian + f"{count}I", blob) if typ == 8: return struct.unpack(endian + f"{count}h", blob) if typ == 12: return struct.unpack(endian + f"{count}d", blob) if typ == 11: return struct.unpack(endian + f"{count}f", blob) raise ValueError(f"unsupported TIFF type {typ}") def read_ifd(path: Path) -> dict[int, tuple]: with path.open("rb") as handle: header = handle.read(8) if len(header) < 8 or header[:2] not in (b"II", b"MM"): raise ValueError(f"{path.name}: not a TIFF (header {header[:4]!r})") endian = "<" if header[:2] == b"II" else ">" magic, ifd_offset = struct.unpack(endian + "HI", header[2:8]) if magic != 42: raise ValueError(f"{path.name}: TIFF version {magic} is not 42") handle.seek(ifd_offset) count_bytes = handle.read(2) if len(count_bytes) != 2: raise ValueError(f"{path.name}: IFD at {ifd_offset} is unreadable") (count,) = struct.unpack(endian + "H", count_bytes) entries: list[tuple[int, int, int, bytes]] = [] for _ in range(count): entry = handle.read(12) if len(entry) != 12: raise ValueError(f"{path.name}: truncated IFD entry") tag, typ, field_count = struct.unpack(endian + "HHI", entry[:8]) entries.append((tag, typ, field_count, entry[8:12])) needed = {TAG_WIDTH, TAG_HEIGHT, TAG_PIXEL_SCALE, TAG_TIEPOINT, TAG_GEOKEY} values: dict[int, tuple] = {} for tag, typ, field_count, raw in entries: if tag not in needed: continue values[tag] = read_values(handle, endian, typ, field_count, raw) return values def projected_epsg(path: Path, geokeys: tuple) -> int: if len(geokeys) < 4: raise ValueError(f"{path.name}: GeoKeyDirectory is shorter than its header") key_count = geokeys[3] needed = 4 + key_count * 4 if len(geokeys) < needed: raise ValueError( f"{path.name}: GeoKeyDirectory has {len(geokeys)} shorts, needs {needed}" ) for index in range(key_count): base = 4 + index * 4 key_id, location, _count, value = geokeys[base : base + 4] if key_id != GEOKEY_PROJECTED_CS: continue if location != 0: raise ValueError( f"{path.name}: ProjectedCSTypeGeoKey location {location} is not an inline value" ) return int(value) raise ValueError(f"{path.name}: missing ProjectedCSTypeGeoKey {GEOKEY_PROJECTED_CS}") def read_blatt(path: Path) -> tuple[Blatt | None, list[str]]: """Return the Blatt, or validation errors that name the file and the bad value.""" errors: list[str] = [] try: tags = read_ifd(path) except (OSError, ValueError, struct.error) as exc: return None, [str(exc) if str(exc).startswith(path.name) else f"{path.name}: {exc}"] def require(tag: int, label: str) -> tuple | None: if tag not in tags: errors.append(f"{path.name}: missing GeoTIFF tag {tag} ({label})") return None return tags[tag] width_values = require(TAG_WIDTH, "ImageWidth") height_values = require(TAG_HEIGHT, "ImageLength") scale_values = require(TAG_PIXEL_SCALE, "ModelPixelScale") tie_values = require(TAG_TIEPOINT, "ModelTiepoint") key_values = require(TAG_GEOKEY, "GeoKeyDirectory") if errors: return None, errors assert width_values and height_values and scale_values and tie_values and key_values width = int(width_values[0]) height = int(height_values[0]) if width != EXPECTED_PIXELS: errors.append(f"{path.name}: width {width} is not {EXPECTED_PIXELS}") if height != EXPECTED_PIXELS: errors.append(f"{path.name}: height {height} is not {EXPECTED_PIXELS}") if len(scale_values) < 2: errors.append(f"{path.name}: pixel scale {scale_values} is not 1.25") scale_x = scale_y = None else: scale_x = float(scale_values[0]) scale_y = float(scale_values[1]) if abs(scale_x - NATIVE_RESOLUTION) > SCALE_TOLERANCE or abs(scale_y - NATIVE_RESOLUTION) > SCALE_TOLERANCE: errors.append(f"{path.name}: pixel scale ({scale_x}, {scale_y}) is not 1.25") if len(tie_values) < 6: errors.append(f"{path.name}: tiepoint {tie_values} does not have six values") return None, errors pixel_i, pixel_j, pixel_k = (float(tie_values[0]), float(tie_values[1]), float(tie_values[2])) if ( abs(pixel_i) > SCALE_TOLERANCE or abs(pixel_j) > SCALE_TOLERANCE or abs(pixel_k) > SCALE_TOLERANCE ): errors.append( f"{path.name}: tiepoint pixel ({pixel_i}, {pixel_j}, {pixel_k}) is not (0, 0, 0)" ) tie_x = float(tie_values[3]) tie_y = float(tie_values[4]) try: epsg = projected_epsg(path, key_values) except ValueError as exc: errors.append(str(exc)) epsg = None if epsg is not None and epsg != EXPECTED_EPSG: errors.append(f"{path.name}: EPSG {epsg} is not {EXPECTED_EPSG}") if scale_x is None or scale_y is None: return None, errors span_e = width * scale_x span_n = height * scale_y if abs(span_e - EXPECTED_SPAN_M) > SPAN_TOLERANCE_M or abs(span_n - EXPECTED_SPAN_M) > SPAN_TOLERANCE_M: errors.append( f"{path.name}: ground span {span_e} m by {span_n} m is not " f"{EXPECTED_SPAN_M} m by {EXPECTED_SPAN_M} m" ) world = path.with_suffix(".tfw") if not world.is_file(): errors.append(f"{path.name}: missing world file {world.name}") if errors: return None, errors # Tiepoint (0, 0) is the upper-left corner. The world file's origin is # the center of that pixel; the tags are enough once they validate. return Blatt( path=path, min_e=tie_x, max_n=tie_y, max_e=tie_x + span_e, min_n=tie_y - span_n, ), [] def load_blaetter(directory: Path) -> tuple[list[Blatt], list[str]]: if not directory.is_dir(): return [], [f"{directory} is not a directory"] paths = sorted(directory.glob("DTK25_*_col.tif")) if not paths: return [], [f"no Blätter matching DTK25_*_col.tif in {directory}"] sheets: list[Blatt] = [] errors: list[str] = [] for path in paths: sheet, sheet_errors = read_blatt(path) errors.extend(sheet_errors) if sheet is not None: sheets.append(sheet) sheets.sort(key=lambda sheet: (sheet.min_n, sheet.min_e, sheet.path.name)) return sheets, errors def write_uebersicht(path: Path, sheets: list[Blatt]) -> None: document = { "tileSizePx": TILE_SIZE_PX, "resolutions": RESOLUTIONS, "blaetter": [sheet.as_json() for sheet in sheets], "gebiete": [{"name": "Hessen", "covers": "all"}], } path.parent.mkdir(parents=True, exist_ok=True) path.write_text(json.dumps(document, indent=2) + "\n", encoding="utf-8") def gdal_available() -> bool: return shutil.which("gdalwarp") is not None and shutil.which("gdal_translate") is not None def write_kachel(incoming, out: Path) -> None: """Store one Kachel. A shared edge keeps opaque pixels from both Blätter. A tile often crosses a Blatt boundary. Writing the second Blatt used to replace the file, so one side of the Kachel stayed transparent and the map showed the background through that strip. """ from PIL import Image layer = incoming.convert("RGBA") if layer is incoming: layer = incoming.copy() try: if out.is_file() and out.stat().st_size > 0: with Image.open(out) as existing: base = existing.convert("RGBA") if base is existing: base = existing.copy() base.alpha_composite(layer) base.save(out, format="PNG") base.close() else: layer.save(out, format="PNG") finally: layer.close() def run_gdal(command: list[str]) -> None: env = os.environ.copy() env["GDAL_PAM_ENABLED"] = "NO" result = subprocess.run(command, capture_output=True, text=True, env=env) if result.returncode != 0: detail = (result.stderr or result.stdout).strip() raise SystemExit(f"{command[0]} failed ({result.returncode}): {detail}") def tile_with_gdal(sheets: list[Blatt], dest: Path) -> int: written = 0 with tempfile.TemporaryDirectory(prefix="kacheln-") as tmp: for sheet in sheets: for z, resolution in enumerate(RESOLUTIONS): ground = resolution * TILE_SIZE_PX cols = list(tile_indices(sheet.min_e, sheet.max_e, ground)) rows = list(tile_indices(sheet.min_n, sheet.max_n, ground)) if not cols or not rows: continue te_min_e = cols[0] * ground te_max_e = (cols[-1] + 1) * ground te_min_n = rows[0] * ground te_max_n = (rows[-1] + 1) * ground warped = Path(tmp) / f"{sheet.path.stem}-z{z}.tif" resample = "near" if z == 0 else "lanczos" run_gdal( [ "gdalwarp", "-q", "-t_srs", f"EPSG:{EXPECTED_EPSG}", "-te", str(te_min_e), str(te_min_n), str(te_max_e), str(te_max_n), "-tr", str(resolution), str(resolution), "-r", resample, "-dstalpha", "-of", "GTiff", "-co", "TILED=YES", "-co", "COMPRESS=DEFLATE", str(sheet.path), str(warped), ] ) tags = read_ifd(warped) width = int(tags[TAG_WIDTH][0]) height = int(tags[TAG_HEIGHT][0]) expect_w = len(cols) * TILE_SIZE_PX expect_h = len(rows) * TILE_SIZE_PX if width != expect_w or height != expect_h: raise SystemExit( f"{sheet.path.name}: gdalwarp size {width} by {height} " f"is not {expect_w} by {expect_h} at {resolution} m" ) for col in cols: column_dir = dest / "kacheln" / str(z) / str(col) column_dir.mkdir(parents=True, exist_ok=True) x_off = int(round((col * ground - te_min_e) / resolution)) for row in rows: y_off = int(round((te_max_n - (row + 1) * ground) / resolution)) out = column_dir / f"{row}.png" fresh = Path(tmp) / f"tile-{z}-{col}-{row}.png" run_gdal( [ "gdal_translate", "-q", "-srcwin", str(x_off), str(y_off), str(TILE_SIZE_PX), str(TILE_SIZE_PX), "-of", "PNG", str(warped), str(fresh), ] ) from PIL import Image with Image.open(fresh) as incoming: write_kachel(incoming, out) fresh.unlink(missing_ok=True) written += 1 print( f" z={z} {resolution} m {len(cols) * len(rows)} tiles {sheet.path.name}", flush=True, ) return written def paste_tile(level, sheet: Blatt, resolution: float, col: int, row: int): from PIL import Image ground = resolution * TILE_SIZE_PX tile_min_e = col * ground tile_max_n = (row + 1) * ground what = f"{sheet.path.name}: tile {resolution} m {col}/{row}" src_x = snap((tile_min_e - sheet.min_e) / resolution, what) src_y = snap((sheet.max_n - tile_max_n) / resolution, what) src_w, src_h = level.size left = max(src_x, 0) upper = max(src_y, 0) right = min(src_x + TILE_SIZE_PX, src_w) lower = min(src_y + TILE_SIZE_PX, src_h) if right <= left or lower <= upper: raise SystemExit(f"{what} intersects the Blatt but has no pixels") tile = Image.new("RGBA", (TILE_SIZE_PX, TILE_SIZE_PX), (0, 0, 0, 0)) crop = level.crop((left, upper, right, lower)) tile.paste(crop, (left - src_x, upper - src_y)) return tile def tile_with_pillow(sheets: list[Blatt], dest: Path) -> int | None: try: from PIL import Image except ImportError as exc: print( "Pillow could not be imported, and gdalwarp/gdal_translate are not on PATH. " f"The Übersicht was written. No tiles were written. ({exc})", file=sys.stderr, ) return None # A Blatt is 9600² pixels, just over Pillow's default bomb limit. Image.MAX_IMAGE_PIXELS = EXPECTED_PIXELS * EXPECTED_PIXELS written = 0 for sheet in sheets: try: with Image.open(sheet.path) as opened: opened.load() native = opened.convert("RGBA") except Exception as exc: print( f"Pillow could not decode {sheet.path.name}: {exc}. " "The Übersicht was written. No tiles were written.", file=sys.stderr, ) return None if native.size != (EXPECTED_PIXELS, EXPECTED_PIXELS): print( f"{sheet.path.name}: Pillow size {native.size[0]} by {native.size[1]} " f"is not {EXPECTED_PIXELS} by {EXPECTED_PIXELS}. " "The Übersicht was written. No tiles were written.", file=sys.stderr, ) native.close() return None try: for z, resolution in enumerate(RESOLUTIONS): scale = NATIVE_RESOLUTION / resolution level_size = ( snap(native.width * scale, f"{sheet.path.name}: resized width at {resolution} m"), snap(native.height * scale, f"{sheet.path.name}: resized height at {resolution} m"), ) level = native if z == 0 else native.resize(level_size, Image.Resampling.LANCZOS) ground = resolution * TILE_SIZE_PX cols = list(tile_indices(sheet.min_e, sheet.max_e, ground)) rows = list(tile_indices(sheet.min_n, sheet.max_n, ground)) level_count = 0 try: for col in cols: column_dir = dest / "kacheln" / str(z) / str(col) column_dir.mkdir(parents=True, exist_ok=True) for row in rows: tile_min_e = col * ground tile_min_n = row * ground if not intersects( tile_min_e, tile_min_n, tile_min_e + ground, tile_min_n + ground, sheet, ): continue out = column_dir / f"{row}.png" tile = paste_tile(level, sheet, resolution, col, row) write_kachel(tile, out) tile.close() written += 1 level_count += 1 finally: if level is not native: level.close() print( f" z={z} {resolution} m {level_count} tiles {sheet.path.name}", flush=True, ) finally: native.close() return written def select_sheets(sheets: list[Blatt], only: str | None, tile_all: bool) -> tuple[list[Blatt], str | None]: if tile_all and only: return [], "pass either --only or --all, not both" if tile_all: return sheets, None if only: name = Path(only).name chosen = [sheet for sheet in sheets if sheet.path.name == name] if not chosen: return [], f"no Blatt named {name}" return chosen, None return [sheets[0]], None def parse_args(argv: list[str] | None) -> argparse.Namespace: parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--karten", type=Path, required=True, help="directory of DTK25 GeoTIFFs") parser.add_argument("--uebersicht", type=Path, required=True, help="Übersicht JSON to write") parser.add_argument( "--tiles", type=Path, help="write {DIR}/kacheln/{z}/{col}/{row}.png; pass public (one Blatt unless --all)", ) parser.add_argument( "--only", help="filename of the single Blatt to tile, e.g. DTK25_464000_5576000_col.tif", ) parser.add_argument("--all", action="store_true", help="tile every Blatt") return parser.parse_args(argv) def main(argv: list[str] | None = None) -> int: args = parse_args(argv) sheets, errors = load_blaetter(args.karten) if errors: for error in errors: print(error, file=sys.stderr) return 1 write_uebersicht(args.uebersicht, sheets) print(f"{len(sheets)} Blätter") print(f"wrote {args.uebersicht}") if args.tiles is None: return 0 chosen, error = select_sheets(sheets, args.only, args.all) if error: print(error, file=sys.stderr) return 1 if not args.all and args.only is None: print(f"tiling one Blatt ({chosen[0].path.name})") if gdal_available(): print("tiling with gdalwarp and gdal_translate") written = tile_with_gdal(chosen, args.tiles) else: print("gdalwarp/gdal_translate not on PATH; tiling with Pillow") written = tile_with_pillow(chosen, args.tiles) if written is None: return 1 print(f"{written} PNG tiles written to {args.tiles / 'kacheln'}") return 0 if __name__ == "__main__": sys.exit(main())