Files
GeoData/export_ortho_tiles.py

111 lines
3.9 KiB
Python

#!/usr/bin/env python3
"""Export orthophoto tiles aligned to the terrain grid."""
from __future__ import annotations
import argparse
import csv
import glob
import os
from osgeo import gdal
from gdal_utils import build_vrt, ensure_dir, ensure_parent, open_dataset, resolve_first_existing
gdal.UseExceptions()
def parse_args() -> argparse.Namespace:
parser = argparse.ArgumentParser(description="Export orthophoto tiles aligned to the terrain grid.")
parser.add_argument("--raw-ortho-dir", default="raw_dop/jp2", help="Legacy directory containing JP2 orthophoto tiles.")
parser.add_argument(
"--raw-ortho-dir-new",
dest="raw_ortho_dir_new",
default="raw/dop20/jp2",
help="Preferred directory containing JP2 orthophoto tiles (new layout).",
)
parser.add_argument(
"--raw-ortho-dir-alt",
dest="raw_ortho_dir_alt",
default="raw/dop/jp2",
help="Alternate directory containing JP2 orthophoto tiles (intermediate layout).",
)
parser.add_argument("--vrt-path", default="work/dop.vrt", help="Path to build/read the orthophoto VRT.")
parser.add_argument("--tile-index", default="export_unity/tile_index.csv", help="Tile manifest from heightmap export.")
parser.add_argument("--out-dir", default="export_unity/ortho_jpg", help="Output directory for cropped orthophotos.")
parser.add_argument(
"--out-res",
type=int,
default=2048,
help="Output resolution per tile (default matches 1000 m tiles at ~0.5 m/px).",
)
parser.add_argument("--jpeg-quality", type=int, default=90, help="JPEG quality for exported tiles.")
return parser.parse_args()
def export_orthos(args: argparse.Namespace) -> int:
ensure_dir("work")
ensure_dir(args.out_dir)
ensure_parent(args.vrt_path)
candidates = [args.raw_ortho_dir_new, args.raw_ortho_dir_alt, args.raw_ortho_dir]
raw_ortho_dir = resolve_first_existing(candidates, "orthophoto input directory")
jp2_paths = sorted(glob.glob(os.path.join(raw_ortho_dir, "*.jp2")))
if not jp2_paths:
raise SystemExit(
f"No JP2 files found in {raw_ortho_dir}. Run scripts/dlscript_dop20.sh (or use the fallback layouts) first."
)
build_vrt(args.vrt_path, jp2_paths)
vrt_ds = open_dataset(args.vrt_path, f"Could not open VRT at {args.vrt_path}")
if not os.path.exists(args.tile_index):
raise SystemExit(f"Tile index missing: {args.tile_index}. Run export_heightmaps.py first.")
with open(args.tile_index, newline="", encoding="utf-8") as f:
reader = csv.DictReader(f)
written = 0
skipped = 0
for row in reader:
try:
tile_id = row["tile_id"]
xmin = float(row["xmin"])
ymin = float(row["ymin"])
xmax = float(row["xmax"])
ymax = float(row["ymax"])
except (KeyError, ValueError) as exc:
print(f"Skipping malformed row {row}: {exc}")
skipped += 1
continue
out_path = os.path.join(args.out_dir, f"{tile_id}.jpg")
opts = gdal.TranslateOptions(
format="JPEG",
width=args.out_res,
height=args.out_res,
projWin=(xmin, ymax, xmax, ymin), # xmin,xmax,ymax,ymin (upper-left origin)
creationOptions=[f"QUALITY={args.jpeg_quality}", "WORLDFILE=YES"],
)
try:
gdal.Translate(out_path, vrt_ds, options=opts)
except RuntimeError as exc:
print(f"Translate failed for {tile_id}: {exc}")
skipped += 1
continue
written += 1
print(f"Wrote {out_path}")
print(f"Summary: wrote {written} orthophoto tiles; skipped {skipped}.")
return 1 if skipped else 0
def main() -> None:
args = parse_args()
raise SystemExit(export_orthos(args))
if __name__ == "__main__":
main()