เมื่อวิศวกรสำรวจคำนวณเรื่องน้ำ – วิเคราะห์ระบบระบายน้ำบริเวณ RESA

#!/usr/bin/env python3
"""
landxml_to_geotiff.py
Convert LandXML TIN surface to GeoTIFF raster (elevation)

Usage:
  python landxml_to_geotiff.py input.xml output.tif [resolution] [epsg] [crop_area.dxf]

Arguments:
  input.xml    : LandXML file containing TIN surface
  output.tif   : Output GeoTIFF file
  resolution   : Grid resolution in metres (default: 0.10)
  epsg         : CRS EPSG code (default: 32647 = UTM Zone 47N Thailand)

Example:
  python landxml_to_geotiff.py NorthRESA.xml NorthRESA_DEM.tif 0.10 32647

Requirements:
  pip install lxml numpy scipy rasterio

Author: based on workflow by Prajuab Riabroy (priabroy.name)
"""

import sys
import numpy as np
from lxml import etree
from scipy.interpolate import LinearNDInterpolator
import rasterio
from rasterio.transform import from_bounds
from rasterio.crs import CRS

# ─── Parse arguments ────────────────────────────────────────────
def usage():
    print(__doc__)
    sys.exit(1)

if len(sys.argv) < 3:
    usage()

input_xml  = sys.argv[1]
output_tif = sys.argv[2]
resolution = float(sys.argv[3]) if len(sys.argv) > 3 else 0.10
epsg       = int(sys.argv[4])   if len(sys.argv) > 4 else 32647

print(f"Input  : {input_xml}")
print(f"Output : {output_tif}")
print(f"Resolution: {resolution} m")
print(f"CRS    : EPSG:{epsg}")

# ─── Parse LandXML ──────────────────────────────────────────────
print("\nParsing LandXML...")
tree = etree.parse(input_xml)
root = tree.getroot()

# Handle namespace — LandXML may or may not have one
ns = root.nsmap.get(None, '')
if ns:
    ns_prefix = f'{{{ns}}}'
else:
    ns_prefix = ''

def tag(name):
    return f'{ns_prefix}{name}'

# Find first Surface
surface = root.find(f'.//{tag("Surface")}')
if surface is None:
    print("ERROR: No <Surface> element found in LandXML file.")
    sys.exit(1)

surface_name = surface.get('name', 'Unknown')
print(f"Surface: {surface_name}")

# ─── Extract Points ──────────────────────────────────────────────
print("Reading points...")
pnts_elem = surface.find(f'.//{tag("Pnts")}')
if pnts_elem is None:
    print("ERROR: No <Pnts> element found.")
    sys.exit(1)

points = {}
for p in pnts_elem.findall(tag('P')):
    pid = int(p.get('id'))
    coords = p.text.strip().split()
    # LandXML point order: N E Z (Northing, Easting, Elevation)
    n, e, z = float(coords[0]), float(coords[1]), float(coords[2])
    points[pid] = (e, n, z)  # store as (X=Easting, Y=Northing, Z)

print(f"  Points: {len(points):,}")

# ─── Extract Faces (Triangles) ──────────────────────────────────
print("Reading faces...")
faces_elem = surface.find(f'.//{tag("Faces")}')
if faces_elem is None:
    print("ERROR: No <Faces> element found.")
    sys.exit(1)

faces = []
for f in faces_elem.findall(tag('F')):
    ids = f.text.strip().split()
    faces.append((int(ids[0]), int(ids[1]), int(ids[2])))

print(f"  Faces : {len(faces):,}")

# ─── Build point arrays ──────────────────────────────────────────
print("\nBuilding interpolator...")
xs = np.array([points[pid][0] for pid in points])
ys = np.array([points[pid][1] for pid in points])
zs = np.array([points[pid][2] for pid in points])

x_min, x_max = xs.min(), xs.max()
y_min, y_max = ys.min(), ys.max()

print(f"  X range: {x_min:.2f} → {x_max:.2f}  ({x_max-x_min:.1f} m)")
print(f"  Y range: {y_min:.2f} → {y_max:.2f}  ({y_max-y_min:.1f} m)")
print(f"  Z range: {zs.min():.3f} → {zs.max():.3f} m")

# ─── Create output grid ──────────────────────────────────────────
cols = int(np.ceil((x_max - x_min) / resolution)) + 1
rows = int(np.ceil((y_max - y_min) / resolution)) + 1
print(f"\nOutput grid: {cols} × {rows} cells at {resolution}m resolution")
print(f"  Estimated size: {cols * rows / 1e6:.1f} million cells")

# Grid coordinates (centre of each cell)
xi = np.linspace(x_min, x_max, cols)
yi = np.linspace(y_max, y_min, rows)  # top-to-bottom (raster convention)
xi_grid, yi_grid = np.meshgrid(xi, yi)

# ─── Interpolate ─────────────────────────────────────────────────
print("\nInterpolating (LinearNDInterpolator)...")
print("  This may take a moment for large surfaces...")

interp = LinearNDInterpolator(list(zip(xs, ys)), zs)
zi = interp(xi_grid, yi_grid)

# Replace NaN (outside TIN boundary) with nodata
nodata = -9999.0
zi[np.isnan(zi)] = nodata

valid_count = np.sum(zi != nodata)
print(f"  Valid cells: {valid_count:,} / {cols*rows:,} ({100*valid_count/(cols*rows):.1f}%)")

# ─── Write GeoTIFF ───────────────────────────────────────────────
print(f"\nWriting GeoTIFF: {output_tif}")
transform = from_bounds(x_min, y_min, x_max, y_max, cols, rows)
crs = CRS.from_epsg(epsg)

with rasterio.open(
    output_tif, 'w',
    driver   = 'GTiff',
    height   = rows,
    width    = cols,
    count    = 1,
    dtype    = np.float32,
    crs      = crs,
    transform= transform,
    nodata   = nodata,
    compress = 'lzw',
) as dst:
    dst.write(zi.astype(np.float32), 1)

print(f"\n✓ Done! Output: {output_tif}")
print(f"  Grid size   : {cols} × {rows}")
print(f"  Resolution  : {resolution} m")
print(f"  CRS         : EPSG:{epsg}")
print(f"  Z min/max   : {zs.min():.3f} / {zs.max():.3f} m")
print(f"\nReady to open in QGIS or use with dem_viewer.html")