#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
LISCAD "SEEDB" XML export: DrukRef03 TM -> DrukRef23 LCC transformation using NTv2.

IMPORTANT - schema note
------------------------
The sample file (testTM.xml) is NOT LandXML. Its root element is <SEEDB>,
and it uses LISCAD's own native export schema:

  <Points><Point coords="E N" elevation="Z" num="1" .../></Points>
  <Lines><StraightLine num="1" pnts="1 2" length="..." ...>
      <PT pos="1">1</PT><PT pos="2">2</PT>
  </StraightLine></Lines>
  <Polygons><Polygon num="10" lines="1 2 3 4 5 6 7 8 9"
                     coords="E N" area="..." .../></Polygons>

Key differences from the original (LandXML-oriented) script, which is why it
transformed 0 coordinates on this file:

  1. There is no <CgPoint>/<Start>/<End>/<Center> text-based coordinate
     anywhere - those tags don't exist in this schema at all.
  2. Coordinates live in the coords="..." ATTRIBUTE of <Point> (and
     <Polygon>), not in element text.
  3. The coordinate order in coords="..." is Easting Northing (X Y), not
     Northing Easting as LandXML/LISCAD's CgPoint convention would be.
  4. Elevation is a separate elevation="..." attribute on <Point>, not a
     third value appended to the coordinate string.
  5. <StraightLine> does not carry inline Start/End coordinates - it
     references point numbers via pnts="p1 p2". <PT> child elements repeat
     the same point numbers for display/running-distance bookkeeping.
  6. <Polygon> area/centroid must be rebuilt by walking its lines="..."
     attribute (a list of StraightLine numbers) back to points via pnts=.

This script transforms every <Point>, then recomputes <StraightLine length=">
and <Polygon area="/coords="> from the transformed point coordinates.

Usage (OSGeo4W Shell):
    python landxml_to_landxml_ntv2_drukref23_LCC.py ^
        "input.xml" ^
        "output.xml" ^
        "D:\\transoformed23\\d03tod23.gsb" ^
        EPSG:11341 ^
        --src-epsg EPSG:5266
"""

import os
import sys
import math
import argparse
import re
import xml.etree.ElementTree as ET

from osgeo import osr, gdal


def detect_xml_encoding(input_xml):
    """Read the XML declaration encoding; default to UTF-8."""
    with open(input_xml, "rb") as f:
        head = f.read(512)
    m = re.search(br'encoding=["\']([^"\']+)["\']', head, re.I)
    return m.group(1).decode("ascii") if m else "utf-8"


def make_srs(epsg):
    srs = osr.SpatialReference()
    err = srs.ImportFromEPSG(int(str(epsg).split(":")[-1]))
    if err != 0:
        raise RuntimeError("Could not load CRS: {}".format(epsg))

    if hasattr(osr, "OAMS_TRADITIONAL_GIS_ORDER"):
        srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)

    return srs


def make_transformation(src_epsg, dst_epsg, grid_file):
    src = make_srs(src_epsg)
    dst = make_srs(dst_epsg)

    # Source DrukRef03 / Bhutan National Grid TM.
    # Confirmed against the <Projection>/<TransverseMercator>/<Ellipsoid>
    # block in testTM.xml: centralMeridian=90, falseEast=250000,
    # falseNorth=0, scaleFactor=1, ellipsoid a=6378137 / invF=298.257222
    # (GRS80).
    src_proj = (
        "+proj=tmerc +lat_0=0 +lon_0=90 +k=1 "
        "+x_0=250000 +y_0=0 +ellps=GRS80 +units=m +no_defs"
    )

    dst_proj = dst.ExportToProj4()
    if not dst_proj:
        raise RuntimeError(
            "GDAL/PROJ could not export {} to PROJ syntax.".format(dst_epsg)
        )
    dst_proj = dst_proj.replace("+type=crs", "").strip()

    pipeline = (
        "+proj=pipeline "
        "+step +inv " + src_proj + " "
        "+step +proj=hgridshift +grids=" + os.path.abspath(grid_file).replace("\\", "/") + " "
        "+step " + dst_proj
    )

    options = osr.CoordinateTransformationOptions()

    if not options.SetOperation(pipeline):
        raise RuntimeError("Could not set the coordinate transformation pipeline.")

    if hasattr(options, "SetBallparkAllowed"):
        options.SetBallparkAllowed(False)

    if hasattr(options, "SetOnlyBest"):
        options.SetOnlyBest(True)

    ct = osr.CreateCoordinateTransformation(src, dst, options)
    if ct is None:
        raise RuntimeError("Could not create coordinate transformation.")

    return ct, pipeline


def parse_en(text):
    """coords="..." attribute is Easting Northing (X Y)."""
    vals = text.strip().split()
    if len(vals) < 2:
        raise ValueError("coords attribute has fewer than 2 values: {!r}".format(text))
    e = float(vals[0])
    n = float(vals[1])
    return e, n


def transform_points(root, ct, stats):
    """
    Transform every <Point coords="E N" elevation="Z" num="..."/>.
    Returns dict: point num (str) -> (new_e, new_n, new_z)
    """
    points_by_num = {}

    for point in root.iter("Point"):
        num = point.get("num")
        coords = point.get("coords")
        if not coords:
            continue

        try:
            e, n = parse_en(coords)
            z_txt = point.get("elevation")
            z = float(z_txt) if z_txt not in (None, "") else 0.0

            x2, y2, z2 = ct.TransformPoint(e, n, z)

            point.set("coords", "{:.6f} {:.6f}".format(x2, y2))
            if z_txt not in (None, ""):
                point.set("elevation", "{:.6f}".format(z2))

            points_by_num[num] = (x2, y2, z2)
            stats["coordinates"] += 1

        except Exception as exc:
            stats["warnings"].append(
                "Could not transform Point num={}: {}".format(num, exc)
            )

    return points_by_num


def update_straightline_lengths(root, points_by_num, stats):
    """Recalculate StraightLine/@length from the transformed endpoints."""
    lines_by_num = {}
    count = 0

    for line in root.iter("StraightLine"):
        num = line.get("num")
        lines_by_num[num] = line

        pnts = (line.get("pnts") or "").split()
        if len(pnts) != 2:
            stats["warnings"].append(
                "StraightLine num={} does not have exactly 2 pnts; length left unchanged".format(num)
            )
            continue

        p1, p2 = pnts
        if p1 not in points_by_num or p2 not in points_by_num:
            stats["warnings"].append(
                "StraightLine num={} references missing point(s) {}/{}".format(num, p1, p2)
            )
            continue

        e1, n1, _ = points_by_num[p1]
        e2, n2, _ = points_by_num[p2]

        line.set("length", "{:.6f}".format(math.hypot(e2 - e1, n2 - n1)))
        count += 1

    return lines_by_num, count


def polygon_area_and_centroid(ordered_en):
    """Shoelace area (E,N order) + simple vertex-average centroid."""
    n_pts = len(ordered_en)
    if n_pts < 3:
        return 0.0, None

    area2 = 0.0
    for i in range(n_pts):
        e1, n1 = ordered_en[i]
        e2, n2 = ordered_en[(i + 1) % n_pts]
        area2 += e1 * n2 - e2 * n1

    area = abs(area2) / 2.0
    cx = sum(p[0] for p in ordered_en) / n_pts
    cy = sum(p[1] for p in ordered_en) / n_pts
    return area, (cx, cy)


def update_polygon_areas(root, points_by_num, lines_by_num, stats):
    """
    Recalculate Polygon/@area and Polygon/@coords by walking each polygon's
    lines="..." attribute back to StraightLine/@pnts -> Point coordinates.

    Polygons whose lines list references anything other than a StraightLine
    (e.g. an arc/curve num that isn't in lines_by_num) are left unchanged
    and reported as a warning, same caution as the LandXML version.
    """
    updated = 0
    skipped = 0

    for polygon in root.iter("Polygon"):
        num = polygon.get("num")
        line_nums = (polygon.get("lines") or "").split()

        if not line_nums:
            continue

        ordered_pts = []
        ok = True

        for i, lnum in enumerate(line_nums):
            line_el = lines_by_num.get(lnum)
            if line_el is None:
                stats["warnings"].append(
                    "Polygon num={} references line {} which is not a StraightLine; "
                    "area/coords left unchanged".format(num, lnum)
                )
                ok = False
                break

            pnts = (line_el.get("pnts") or "").split()
            if len(pnts) != 2:
                ok = False
                break

            p1, p2 = pnts
            if i == 0:
                ordered_pts.append(p1)
            ordered_pts.append(p2)

        if not ok:
            skipped += 1
            continue

        try:
            en_list = [points_by_num[p][:2] for p in ordered_pts]
        except KeyError as exc:
            stats["warnings"].append(
                "Polygon num={} references missing point {}".format(num, exc)
            )
            skipped += 1
            continue

        area, centroid = polygon_area_and_centroid(en_list)
        polygon.set("area", "{:.6f}".format(area))
        if centroid is not None:
            polygon.set("coords", "{:.6f} {:.6f}".format(centroid[0], centroid[1]))
        updated += 1

    return updated, skipped


def transform_seedb(input_xml, output_xml, grid_file,
                     src_epsg="EPSG:5266", dst_epsg="EPSG:11341",
                     update_measurements=True):

    if not os.path.isfile(input_xml):
        raise FileNotFoundError("Input XML not found: {}".format(input_xml))

    if not os.path.isfile(grid_file):
        raise FileNotFoundError("NTv2 grid not found: {}".format(grid_file))

    output_dir = os.path.dirname(os.path.abspath(output_xml))
    if output_dir:
        os.makedirs(output_dir, exist_ok=True)

    print("=" * 72)
    print("LISCAD SEEDB XML: DrukRef03 TM -> DrukRef23 LCC")
    print("=" * 72)
    print("Input XML : {}".format(input_xml))
    print("Output XML: {}".format(output_xml))
    print("Grid      : {}".format(grid_file))
    print("Source CRS: {}".format(src_epsg))
    print("Target CRS: {}".format(dst_epsg))
    print()

    print("GDAL version: {}".format(gdal.VersionInfo("--version")))

    ct, pipeline = make_transformation(src_epsg, dst_epsg, grid_file)

    print("NTv2 transformation pipeline:")
    print(pipeline)
    print()

    tree = ET.parse(input_xml)
    root = tree.getroot()

    if root.tag != "SEEDB":
        print("WARNING: root element is <{}>, not <SEEDB> - this script targets "
              "the LISCAD SEEDB export schema; results may be incomplete.".format(root.tag))

    stats = {"coordinates": 0, "warnings": []}

    points_by_num = transform_points(root, ct, stats)

    line_count = 0
    poly_count = 0
    poly_skipped = 0

    if update_measurements:
        lines_by_num, line_count = update_straightline_lengths(root, points_by_num, stats)
        poly_count, poly_skipped = update_polygon_areas(root, points_by_num, lines_by_num, stats)

    output_encoding = detect_xml_encoding(input_xml)
    tree.write(output_xml, encoding=output_encoding, xml_declaration=True)

    print("Transformation completed.")
    print("Points transformed       : {}".format(stats["coordinates"]))
    print("StraightLine lengths set : {}".format(line_count))
    print("Polygon areas set        : {}".format(poly_count))
    if poly_skipped:
        print("Polygon areas skipped    : {}".format(poly_skipped))

    if stats["warnings"]:
        print()
        print("WARNINGS:")
        for w in stats["warnings"]:
            print(" - {}".format(w))

    print()
    print("NOTE: The <Projection>/<TransverseMercator>/<Ellipsoid> metadata block")
    print("      still describes the source DrukRef03 TM definition. This script")
    print("      deliberately leaves it untouched, since rewriting it to an LCC")
    print("      definition uses a different LISCAD projection element and could")
    print("      make the file unreadable by LISCAD. Reset the projection to")
    print("      DrukRef23 LCC / BNG2023 (EPSG:11341) manually inside LISCAD")
    print("      after import, or via Job Properties.")

    print()
    print("Output:")
    print(os.path.abspath(output_xml))
    print("=" * 72)


def main():
    parser = argparse.ArgumentParser(
        description="Transform a LISCAD SEEDB XML export from DrukRef03 TM to DrukRef23 LCC using NTv2."
    )

    parser.add_argument("input_xml", help="Input SEEDB XML file")
    parser.add_argument("output_xml", help="Output SEEDB XML file")
    parser.add_argument("grid_file", help="NTv2 grid file, e.g. d03tod23.gsb")
    parser.add_argument("dst_epsg", help="Target CRS, e.g. EPSG:11341")
    parser.add_argument(
        "--src-epsg",
        default="EPSG:5266",
        help="Source CRS (default: EPSG:5266)"
    )
    parser.add_argument(
        "--no-update-measurements",
        action="store_true",
        help="Do not recalculate StraightLine/@length and Polygon/@area,@coords"
    )

    args = parser.parse_args()

    try:
        transform_seedb(
            args.input_xml,
            args.output_xml,
            args.grid_file,
            args.src_epsg,
            args.dst_epsg,
            update_measurements=not args.no_update_measurements
        )
    except Exception as exc:
        print()
        print("ERROR: {}".format(exc))
        sys.exit(1)


if __name__ == "__main__":
    main()