#! /usr/bin/env python
# filename: leicaxml.py derived from landxml.py
# Copyright 2015 Damien Gaignon <damien.gaignon@gmail.com>
# Copyright 2026 Stefano Costa <steko@iosa.it>
#
# This file is part of Total Open Station.
#
# Total Open Station is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License as
# published by the Free Software Foundation, either version 3 of the
# License, or (at your option) any later version.
#
# Total Open Station is distributed in the hope that it will be
# useful, but WITHOUT ANY WARRANTY; without even the implied warranty
# of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with Total Open Station.  If not, see
# <http://www.gnu.org/licenses/>.

import math
import re
import time
import xml.etree.ElementTree as xml

try:
    from defusedxml import ElementTree as safe_xml
except ImportError:  # pragma: no cover - optional dependency
    safe_xml = xml

from topaze.totalopenstation.formats.conversion import vertical_to_zenithal

from . import UNKNOWN_POINT, Feature, Parser, Point
from .polar import BasePoint, PolarPoint

# Template string
TEMPLATE = """<?xml version="1.0"?>
                <LandXML xmlns="http://www.landxml.org/schema/LandXML-1.2" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.landxml.org/schema/LandXML-1.2 http://www.landxml.org/schema/LandXML-1.2/LandXML-1.2.xsd" date="" time="" version="1.2">
                    <Units>
                        <Metric areaUnit="squareMeter" linearUnit="meter" volumeUnit="cubicMeter" temperatureUnit="celsius" pressureUnit="milliBars" angularUnit="grads" directionUnit="grads"></Metric>
                    </Units>
                    <Project name="Template"></Project>
                    <Application name="TotalOpen Station" desc="TOPS" manufacturer="" version="" manufacturerURL="http://tops.iosa.it/" timeStamp=""></Application>
                </LandXML>"""

DEFAULT_NS = "http://www.landxml.org/schema/LandXML-1.2"
DECLARATION = '<?xml version="1.0" encoding="UTF-8"?>'


def _indent(elem, level=0):
    """
    A function to indent a XML Element for pretty printing
    :param elem: The element to parse
    :param level: The level of the element in the hierarchy
    :return: The element ready to be pretty printed

    Note : Taken from https://stackoverflow.com/questions/749796/pretty-printing-xml-in-python/4590052#4590052
    """
    tab = "\t"
    i = f"\n{(level + 1)*tab}"
    j = f"\n{level*tab}"
    count = 1
    if len(elem):
        if elem.text is None or elem.text.strip() is None:
            elem.text = i
        for subelem in elem:
            _indent(subelem, level + 1)
            if subelem.tail is None or subelem.tail.strip() is None:
                subelem.tail = i
            elif re.match(r"\s", subelem.tail):
                subelem.tail = i
            count += 1
        if count == len(elem) + 1:
            subelem.tail = j
    return elem


class Survey:
    """
    Populate the survey tag of LandXML.

    Tags should be in order (convention):
        - SurveyHeader
        - Equipment
        - CgPoints
        - InstrumentSetup0
        - InstrumentSetupX
        - ObservationGroup0
        - ObservationGroupX

    """

    def __init__(self):
        """
        Initialize the Survey header tag
        """
        self.survey = xml.Element("Survey")
        xml.SubElement(self.survey, "SurveyHeader", name="from TOPS")
        self.id = 0

    def _tag_position(self, tag):
        """
        Search the position of the tag to be in the right order
        """
        mainTags = ("SurveyHeader", "Equipment", "CgPoints", "InstrumentSetup")
        if tag in mainTags:
            pos = 0

            for index in range(mainTags.index(tag)):
                if self.survey.find(f"./{mainTags[index]}") is not None:
                    pos += 1
                else:
                    pos = 1

            if tag == "InstrumentSetup":
                pos += len(self.survey.findall("./InstrumentSetup"))

            return pos
        if tag == "CgPoint":
            cgpoints = self.survey.find("CgPoints")
            pos = len(cgpoints.findall("./CgPoint"))
            return pos

    def equipment(self, **kwargs):
        """
        Populate the Equipment tag.

        No attribut is mandatory so kwargs can be empty.

        Arguments in kwargs for Equipment:
            TODO
        """
        # kwargs = {key: str(value) if value is not None else value for key,value in kwargs.items()}
        # Verification of position of the element
        pos = self._tag_position("Equipment")

        # Creation of Equipment tag, subelement of Survey
        equipment = xml.Element("Equipment")
        self.survey.insert(pos, equipment)

    def cg_point(self, **kwargs):
        """
        Populate the CgPoints and CgPoint tags.

        No attribut is mandatory so kwargs can be empty.

        Arguments in kwargs for CgPoint:
            - point_name -> name          attrib  of CgPoint part of CgPoints
            - pid        -> pntRef        attrib  of CgPoint part of CgPoints
            - x          -> x coordinate  element of CgPoint part of CgPoints
            - y          -> y coordinate  element of CgPoint part of CgPoints
            - z          -> z coordinate  element of CgPoint part of CgPoints
            - attrib     -> attribX       attrib  of Property part of Feature

        All Feature tags should be after Cgpoint ones
        """

        # kwargs = {key: str(value) if value is not None else value for key,value in kwargs.items()}
        # Verification of position of the element
        pos = self._tag_position("CgPoints")

        # Creation of CgPoints tag, subelement of Survey if it does not exist
        cgpoints = self.survey.find("./CgPoints")
        if cgpoints is None:
            cgpoints = xml.Element("CgPoints")
            self.survey.insert(pos, cgpoints)

        # Creation of CgPoint tag, subelement of CgPoints
        cgpoint = xml.Element("CgPoint")
        cgpoints.insert(self._tag_position("CgPoint"), cgpoint)
        # Fill of CgPoint attributes
        if "point_name" in kwargs:
            cgpoint.set("name", str(kwargs["point_name"]))
        if "pid" in kwargs:
            cgpoint.set("pntRef", str(kwargs["pid"]))
        if "x" in kwargs:
            cgpoint.text = f"{kwargs['x']} {kwargs['y']}"
            if "z" in kwargs:
                cgpoint.text += f" {kwargs['z']}"
        # attrib is not mandatory in CgPoints so this is a feature
        if "attrib" in kwargs:
            cgpoint.set("featureRef", f"feature{kwargs['point_name']}")
            feature = xml.Element("Feature")
            feature.set("name", f"feature{kwargs['point_name']}")
            # feature_property
            for i in range(len(kwargs["attrib"])):
                xml.SubElement(
                    feature,
                    "Property",
                    label=f"attrib{i + 1}",
                    value=str(kwargs["attrib"][i]),
                )
            cgpoints.append(feature)

    def setup(self, **kwargs):
        """
        Populate the InstrumentSetup and ObservationGroup tags.

        InstrumentSetup is not a mandatory tag.

        Some attributes are mandatory but can be empty (so kwargs can be empty).

        Arguments in kwargs for InstrumentSetup (attrib mandatory as attrib*):
            - point_name -> stationName      attrib* of InstrumentSetup
            - ih         -> instrumentHeight attrib* of InstrumentSetup
            - pid        -> pntRef           attrib  of InstrumentPoint part of InstrumentSetup
            - instru_x   -> x coordinate     element of InstrumentPoint part of InstrumentSetup
            - instru_y   -> y coordinate     element of InstrumentPoint part of InstrumentSetup
            - instru_z   -> z coordinate     element of InstrumentPoint part of InstrumentSetup
            - attrib     -> attribX          attrib  of Property part of Feature

        Arguments in kwargs for Observationgroup (attrib mandatory as attrib*):
            - circle    -> azimuth          attrib* of Backsight part of ObservationGroup
            - back_name -> name             attrib  of BacksightPoint part of Backsight
            - back_x    -> x coordinate     element of BacksightPoint part of Backsight
            - back_y    -> y coordinate     element of BacksightPoint part of Backsight
            - back_z    -> z coordinate     element of BacksightPoint part of Backsight
        """

        # kwargs = {key: str(value) if value is not None else value for key,value in kwargs.items()}
        # Verification of position of the element
        pos = self._tag_position("InstrumentSetup")

        # Creation of InstrumentSetup tag, subelement of Survey
        instrument_setup = xml.Element(
            "InstrumentSetup", id=f"setup{self.id}", stationName="", instrumentHeight=""
        )
        # Fill of InstrumentSetup attributes
        if "point_name" in kwargs:
            instrument_setup.set("stationName", str(kwargs["point_name"]))
        if "ih" in kwargs:
            instrument_setup.set("instrumentHeight", str(kwargs["ih"]))
        if "hz0" in kwargs:
            instrument_setup.set("orientationAzimuth", str(kwargs["hz0"]))
        # attrib is not mandatory in InstrumentSetup so this is a feature
        if "attrib" in kwargs and kwargs["attrib"]:
            if instrument_setup.find("./Feature") is None:
                feature = xml.SubElement(instrument_setup, "Feature")
            # feature_property
            for i in range(len(kwargs["attrib"])):
                xml.SubElement(
                    feature,
                    "Property",
                    label=f"attrib{i + 1}",
                    value=str(kwargs["attrib"][i]),
                )

        # Creation of InstrumentPoint tag, subelement of InstrumentSetup
        instrument_point = xml.SubElement(instrument_setup, "InstrumentPoint")
        # Fill of InstrumentPoint attributes
        if "pid" in kwargs:
            instrument_point.set("pntRef", str(kwargs["pid"]))
        if "instru_x" in kwargs:
            instrument_point.text = (
                f"{kwargs['instru_x']} {kwargs['instru_y']} {kwargs['instru_z']}"
            )
        # instrument_setup.append(instrument_point)
        self.survey.insert(pos, instrument_setup)

        # Creation of ObservationGroup tag, subelement of Survey
        observation_group = xml.Element(
            "ObservationGroup", id=f"o{self.id}", setupID=f"setup{self.id}"
        )
        # Creation of Backsight tag, subelement of ObservationGroup
        if "circle" in kwargs or "back_x" in kwargs:
            backsight = xml.SubElement(observation_group, "Backsight", circle="0.")
            # Fill of Backsight attributes
            if "circle" in kwargs:
                backsight.set("circle", str(kwargs["circle"]))
            if "back_x" in kwargs:
                # Creation of BacksightPoint tag, subelement of Backsight
                backsight_point = xml.SubElement(backsight, "BacksightPoint")
                # Fill of BacksightPoint attributes
                if "back_name" in kwargs:
                    backsight_point.set("name", str(kwargs["back_name"]))
                if "back_x" in kwargs:
                    backsight_point.text = (
                        f"{kwargs['back_x']} {kwargs['back_y']} {kwargs['back_z']}"
                    )
        self.survey.insert(pos * 2, observation_group)

        # ID can be raise
        self.id += 1

    def raw_observation(self, **kwargs):
        """
        Populate the RawObservation tag.

        No attribut is mandatory so kwargs can be empty.

        Arguments in kwargs for RawObservation:
            - th              -> targetHeight        attrib  of RawObservation
            - angle           -> horizAngle          attrib  of RawObservation
            - z_angle         -> zenithAngle         attrib  of RawObservation
            - dist            -> Distance            attrib  of RawObservation
            - dist_type       -> Type of distance    attrib  of RawObservation
            - point_name      -> desc                attrib  of TargetPoint part of RawObservation
            - pid             -> pntRef              attrib  of TargetPoint part of RawObservation
            - x               -> x coordinate        element of TargetPoint part of RawObservation
            - y               -> y coordinate        element of TargetPoint part of RawObservation
            - z               -> z coordinate        element of TargetPoint part of RawObservation
            - ih              -> instrumentHeight    attrib  of Property part of Feature
            - ppm             -> edmAccuracyppm      attrib  of Property part of Feature
            - prism_constant  -> edmAccuracyConstant attrib  of Property part of Feature
            - attrib          -> attribX             attrib  of Property part of Feature
        """

        # kwargs = {key: str(value) if value is not None else value for key,value in kwargs.items()}
        # When creating a RawObservation tag, it should verified that an ObservationGroup tag exists
        if self.survey.find("./ObservationGroup[@id='o0']") is None:
            self.setup()
        observation_group = self.survey.find(
            f"./ObservationGroup[@id='o{self.id - 1}']"
        )
        # Creation of RawObservation tag, subelement of ObservationGroup
        raw_observation = xml.SubElement(observation_group, "RawObservation")
        # Fill of RawObservation attributes
        if "th" in kwargs:
            raw_observation.set("targetHeight", str(kwargs["th"]))
        if "angle" in kwargs and kwargs["angle"] is not None:
            raw_observation.set("horizAngle", str(kwargs["angle"]))
        if "azimuth" in kwargs and kwargs["azimuth"] is not None:
            raw_observation.set("azimuth", str(kwargs["azimuth"]))
        if "z_angle" in kwargs and kwargs["z_angle"] is not None:
            if kwargs["z_angle_type"] == "dh":
                raw_observation.set("vertDistance", str(kwargs["z_angle"]))
            if kwargs["z_angle_type"] == "z":
                raw_observation.set("zenithAngle", str(kwargs["z_angle"]))
            if kwargs["z_angle_type"] == "v":
                raw_observation.set(
                    "zenithAngle",
                    str(vertical_to_zenithal(kwargs["z_angle"], kwargs["angle_unit"])),
                )
        if "dist" in kwargs and kwargs["dist"] is not None:
            if kwargs["dist_type"] == "s":
                raw_observation.set("slopeDistance", str(kwargs["dist"]))
            if kwargs["dist_type"] == "h":
                raw_observation.set("horizDistance", str(kwargs["dist"]))
        # Creation of TargetPoint tag, subelement of RawObservation
        target_point = xml.SubElement(raw_observation, "TargetPoint")
        # Fill of TargetPoint attributes
        if "point_name" in kwargs:
            target_point.set("desc", str(kwargs["point_name"]))
        if "pid" in kwargs:
            target_point.set("pntRef", str(kwargs["pid"]))
        if "x" in kwargs:
            target_point.text = f"{kwargs['x']} {kwargs['y']}"
            if "z" in kwargs:
                target_point.text += f" {kwargs['z']}"
        # targetHeight is not mandatory in RawObservation so this is a feature
        if "ih" in kwargs and kwargs["ih"] is not None:
            if raw_observation.find("./Feature") is None:
                feature = xml.SubElement(raw_observation, "Feature")
            # feature_property
            xml.SubElement(
                feature, "Property", label="instrumentHeight", value=str(kwargs["ih"])
            )
        # ppm or prism_constant are not mandatory in RawObservation so this is a feature
        if "ppm" in kwargs and kwargs["ppm"] is not None:
            if raw_observation.find("./Feature") is None:
                feature = xml.SubElement(raw_observation, "Feature")
            # feature_property
            xml.SubElement(
                feature, "Property", label="edmAccuracyppm", value=str(kwargs["ppm"])
            )
            xml.SubElement(
                feature,
                "Property",
                label="edmAccuracyConstant",
                value=str(kwargs["prism_constant"]),
            )
        # attrib is not mandatory in RawObservation so this is a feature
        if "attrib" in kwargs:
            if raw_observation.find("./Feature") is None:
                feature = xml.SubElement(raw_observation, "Feature")
            # feature_property
            for i in range(len(kwargs["attrib"])):
                xml.SubElement(
                    feature,
                    "Property",
                    label=f"attrib{i + 1}",
                    value=str(kwargs["attrib"][i]),
                )

    def to_string(self):
        """
        :return:
        """
        return xml.tostring(self.survey)


class LandXML:
    """
    Create the LandXML file.
    """

    def __init__(self):
        xml.register_namespace("", DEFAULT_NS)
        tree = xml.ElementTree(xml.fromstring(TEMPLATE))
        self.root = tree.getroot()

    def append(self, xml_data):
        self.root.append(xml_data)

    def to_string(self):
        cdate = time.strftime("%Y-%m-%d")
        ctime = time.strftime("%H:%M:%S")
        self.root.set("date", cdate)
        self.root.set("time", ctime)
        application = self.root.find(f"{{{DEFAULT_NS}}}Application")
        application.set("timeStamp", f"{cdate}T{ctime}")
        pretty_xml = _indent(self.root)
        return xml.tostring(pretty_xml).decode()


class FormatParser(Parser):
    """Parse survey data stored in a LandXML 1.2 document.

    The parser accepts the two layouts commonly produced by survey software:

    * ``CgPoints`` directly below ``LandXML`` (Leica Captivate, Civil 3D...);
    * ``CgPoints`` below ``Survey`` (historical Total Open Station layout).

    Coordinates in LandXML are always stored as ``northing easting [elevation]``.
    They are exposed to Topaze as ``x=easting, y=northing, z=elevation``. Linear
    values are converted to metres and angular values to grads before creating
    :class:`Feature` objects.
    """

    _LINEAR_FACTORS = {
        "meter": 1.0,
        "metre": 1.0,
        "meters": 1.0,
        "metres": 1.0,
        "millimeter": 0.001,
        "millimetre": 0.001,
        "millimeters": 0.001,
        "millimetres": 0.001,
        "centimeter": 0.01,
        "centimetre": 0.01,
        "centimeters": 0.01,
        "centimetres": 0.01,
        "kilometer": 1000.0,
        "kilometre": 1000.0,
        "kilometers": 1000.0,
        "kilometres": 1000.0,
        "foot": 0.3048,
        "feet": 0.3048,
        "internationalfoot": 0.3048,
        "ussurveyfoot": 1200.0 / 3937.0,
        "ussurveyfeet": 1200.0 / 3937.0,
        "inch": 0.0254,
        "inches": 0.0254,
        "yard": 0.9144,
        "yards": 0.9144,
        "mile": 1609.344,
        "miles": 1609.344,
    }

    _REFERENCE_PURPOSES = {
        "backsight",
        "back sight",
        "resection",
        "orientation",
        "reference",
        "control",
        "control point",
    }

    def __init__(self, data, **options):
        del options
        self.line = safe_xml.fromstring(data)
        if self.line.tag.rsplit("}", 1)[-1] != "LandXML":
            raise ValueError("The document root is not a LandXML element")
        self._namespace = self._namespace_uri(self.line.tag)
        self._linear_unit = "meter"
        self._angle_unit = "radians"
        self._direction_unit = "radians"
        self._linear_factor = 1.0
        self._features_by_name = {}
        for feature in self.line.iter(self._q("Feature")):
            name = feature.attrib.get("name")
            if name:
                self._features_by_name[name] = feature
        self._read_units()

    @staticmethod
    def _namespace_uri(tag):
        if tag.startswith("{") and "}" in tag:
            return tag[1:].split("}", 1)[0]
        return ""

    def _q(self, tag):
        if self._namespace:
            return "{{{}}}{}".format(self._namespace, tag)
        return tag

    @staticmethod
    def _normalise_unit(unit):
        return re.sub(r"[^a-z0-9]", "", (unit or "").lower())

    @staticmethod
    def _format_number(value):
        if value is None:
            return None
        return format(float(value), ".15g")

    def _read_units(self):
        units = self.line.find(self._q("Units"))
        if units is None:
            return

        unit_node = units.find(self._q("Metric"))
        if unit_node is None:
            unit_node = units.find(self._q("Imperial"))
        if unit_node is None:
            return

        self._linear_unit = unit_node.attrib.get("linearUnit", "meter")
        self._angle_unit = unit_node.attrib.get("angularUnit", "radians")
        self._direction_unit = unit_node.attrib.get("directionUnit", self._angle_unit)
        unit_key = self._normalise_unit(self._linear_unit)
        self._linear_factor = self._LINEAR_FACTORS.get(unit_key, 1.0)

    def _linear_value(self, value):
        if value is None or value == "":
            return None
        return float(value) * self._linear_factor

    @staticmethod
    def _packed_dms_to_degrees(value):
        sign = -1.0 if value < 0 else 1.0
        value = abs(value)
        degrees = math.floor(value)
        minutes_seconds = (value - degrees) * 100.0
        minutes = math.floor(minutes_seconds)
        seconds = (minutes_seconds - minutes) * 100.0
        return sign * (degrees + minutes / 60.0 + seconds / 3600.0)

    def _angle_to_gon(self, value, unit=None):
        if value is None or value == "":
            return None

        angle = float(value)
        unit_key = self._normalise_unit(unit or self._angle_unit)
        if unit_key in {"grad", "grads", "gon", "gons", "gradian", "gradians"}:
            result = angle
        elif unit_key in {
            "degree",
            "degrees",
            "decimaldegree",
            "decimaldegrees",
            "deg",
        }:
            result = angle * 10.0 / 9.0
        elif unit_key in {
            "degreesminutesseconds",
            "degreeminutesecond",
            "decimalddmmss",
            "decimaldms",
            "ddmmss",
            "dms",
        }:
            result = self._packed_dms_to_degrees(angle) * 10.0 / 9.0
        elif unit_key in {"mil", "mils"}:
            result = angle / 16.0
        else:
            # LandXML angular and direction values default to radians.
            result = angle * 200.0 / math.pi
        return result

    def _direction_to_gon(self, value):
        return self._angle_to_gon(value, self._direction_unit)

    def _point_from_text(self, text):
        if not text or not text.strip():
            return None
        coords = text.split()
        if len(coords) < 2:
            return None
        try:
            northing = self._linear_value(coords[0])
            easting = self._linear_value(coords[1])
            elevation = self._linear_value(coords[2]) if len(coords) > 2 else None
        except (TypeError, ValueError):
            return None

        if elevation is None:
            return Point(easting, northing)
        return Point(easting, northing, elevation)

    @staticmethod
    def _point_key_candidates(element):
        if element is None:
            return []
        keys = []
        for attr in ("name", "pntRef", "oID", "desc"):
            value = element.attrib.get(attr)
            if value and value not in keys:
                keys.append(value)
        return keys

    def _feature_properties(self, element):
        properties = []
        if element is None:
            return properties
        features = list(element.findall(self._q("Feature")))
        feature_ref = element.attrib.get("featureRef")
        if feature_ref:
            referenced = self._features_by_name.get(feature_ref)
            if referenced is not None and referenced not in features:
                features.append(referenced)
        for feature in features:
            for prop in feature.findall(self._q("Property")):
                label = prop.attrib.get("label", "")
                value = prop.attrib.get("value", "")
                if label or value:
                    properties.append((label, value))
        return properties

    def _terrain_code(self, element):
        if element is None:
            return None
        code = element.attrib.get("code")
        if code is not None and code != "":
            return code

        properties = self._feature_properties(element)
        for label, value in properties:
            if self._normalise_unit(label) in {
                "code",
                "pointcode",
                "featurecode",
                "fieldcode",
            }:
                return value

        # Total Open Station's own LandXML writer historically stores point
        # attributes as ``attrib1``, ``attrib2``... Feature properties.  Keep
        # these values as a space-separated Topaze terrain code when no
        # explicit LandXML ``code`` property is present.
        attribute_values = []
        for label, value in properties:
            label_key = self._normalise_unit(label)
            if value and re.fullmatch(r"(?:attrib|attribute)\d*", label_key):
                attribute_values.append(value)
        return " ".join(attribute_values) or None

    def _cgpoint_containers(self, survey):
        containers = []
        root_cgpoints = self.line.find(self._q("CgPoints"))
        if root_cgpoints is not None:
            containers.append(root_cgpoints)
        if survey is not None:
            survey_cgpoints = survey.find(self._q("CgPoints"))
            if survey_cgpoints is not None and survey_cgpoints not in containers:
                containers.append(survey_cgpoints)
        return containers

    def _build_cgpoint_index(self, survey):
        records = []
        by_name = {}
        by_oid = {}

        for container in self._cgpoint_containers(survey):
            for element in container.findall(self._q("CgPoint")):
                record = {
                    "element": element,
                    "name": element.attrib.get("name") or element.attrib.get("oID"),
                    "oid": element.attrib.get("oID"),
                    "role": (element.attrib.get("role") or "").strip().lower(),
                    "code": self._terrain_code(element),
                    "point": self._point_from_text(element.text),
                }
                records.append(record)
                if record["name"]:
                    by_name.setdefault(record["name"], []).append(record)
                if record["oid"]:
                    by_oid.setdefault(record["oid"], []).append(record)

        def preference(record):
            return (
                record["point"] is not None,
                record["role"] == "control point",
                record["name"] == record["oid"],
            )

        exact = {
            name: max(name_records, key=preference)
            for name, name_records in by_name.items()
        }
        preferred_oid = {
            oid: max(oid_records, key=preference) for oid, oid_records in by_oid.items()
        }
        return records, exact, preferred_oid

    @staticmethod
    def _record_for_key(key, exact, preferred_oid):
        if not key:
            return None
        return exact.get(key) or preferred_oid.get(key)

    def _resolve_record(self, element, exact, preferred_oid, prefer_control=False):
        if element is None:
            return None

        if prefer_control:
            attributes = ("name", "pntRef", "oID", "desc")
        else:
            # For measured targets, pntRef is more precise than name. Leica may
            # use the name of a control point while pntRef identifies a distinct
            # measured occurrence (for example ``S.230@154``).
            attributes = ("pntRef", "name", "oID", "desc")

        candidates = []
        for attribute in attributes:
            key = element.attrib.get(attribute)
            if key and key not in candidates:
                candidates.append(key)

        records = []
        for key in candidates:
            record = self._record_for_key(key, exact, preferred_oid)
            if record is not None and record not in records:
                records.append(record)
                # A measured occurrence may use an alias as ``name`` while its
                # ``oID`` points to the authoritative control point. Follow
                # that relation when reference/control coordinates are wanted.
                if prefer_control and record["oid"]:
                    control_record = preferred_oid.get(record["oid"])
                    if control_record is not None and control_record not in records:
                        records.append(control_record)
        if prefer_control:
            for record in records:
                if record["role"] == "control point" and record["point"] is not None:
                    return record
        for record in records:
            if record["point"] is not None:
                return record
        return records[0] if records else None

    def _setup_id_from_group(self, group):
        setup_id = group.attrib.get("setupID")
        if setup_id:
            return setup_id
        for tag in ("RawObservation", "ReducedObservation", "Backsight"):
            child = group.find(self._q(tag))
            if child is not None and child.attrib.get("setupID"):
                return child.attrib["setupID"]
        group_id = group.attrib.get("id", "")
        if group_id.startswith("o"):
            return group_id[1:]
        return group_id

    def _backsights(self, survey, exact, preferred_oid):
        by_setup = {}
        if survey is None:
            return by_setup

        for group in survey.findall(self._q("ObservationGroup")):
            setup_id = self._setup_id_from_group(group)
            entries = []
            for backsight in group.findall(self._q("Backsight")):
                target = backsight.find(self._q("BacksightPoint"))
                record = self._resolve_record(
                    target, exact, preferred_oid, prefer_control=True
                )
                # The CgPoint contains the authoritative control coordinates.
                # Leica may write a synthetic elevation in BacksightPoint when
                # the reference is only known in planimetry.
                point = record["point"] if record is not None else None
                if point is None and target is not None:
                    point = self._point_from_text(target.text)
                name = None
                if target is not None:
                    name = target.attrib.get("name") or target.attrib.get("desc")
                    if not name:
                        name = target.attrib.get("pntRef")
                entry = {
                    "element": backsight,
                    "target": target,
                    "name": name,
                    "keys": self._point_key_candidates(target),
                    "point": point,
                    "record": record,
                    "used": False,
                }
                if record is not None:
                    for key in (record["name"], record["oid"]):
                        if key and key not in entry["keys"]:
                            entry["keys"].append(key)
                entries.append(entry)
            by_setup.setdefault(setup_id, []).extend(entries)
        return by_setup

    def _matching_backsight(self, setup_id, target, backsights):
        entries = backsights.get(setup_id, [])
        target_keys = self._point_key_candidates(target)
        for entry in entries:
            if set(target_keys) & set(entry["keys"]):
                entry["used"] = True
                return entry
        target_name = target.attrib.get("name") if target is not None else None
        for entry in entries:
            if target_name and entry["name"] == target_name:
                entry["used"] = True
                return entry
        return None

    def _group_station_data(self, groups, exact, preferred_oid):
        for group in groups:
            target = group.find(self._q("TargetPoint"))
            if target is None:
                continue
            name = (
                target.attrib.get("name")
                or target.attrib.get("desc")
                or target.attrib.get("pntRef")
            )
            point = self._point_from_text(target.text)
            record = self._resolve_record(target, exact, preferred_oid)
            if point is None and record is not None:
                point = record["point"]
            code = record["code"] if record is not None else None
            if code is None:
                code = self._terrain_code(target)
            if name or point is not None:
                return {"name": name, "point": point, "code": code}
        return {"name": None, "point": None, "code": None}

    def _station_feature(self, setup, groups, exact, preferred_oid, pid):
        setup_id = setup.attrib.get("id", "")
        instrument_point = setup.find(self._q("InstrumentPoint"))
        group_data = self._group_station_data(groups, exact, preferred_oid)
        station_name = setup.attrib.get("stationName")
        if not station_name and instrument_point is not None:
            station_name = instrument_point.attrib.get(
                "name"
            ) or instrument_point.attrib.get("pntRef")
        if not station_name:
            station_name = group_data["name"]
        if not station_name:
            station_name = setup_id or "station_{}".format(pid)

        point = self._point_from_text(
            instrument_point.text if instrument_point is not None else None
        )
        record = self._resolve_record(instrument_point, exact, preferred_oid)
        if point is None and record is not None:
            point = record["point"]
        if point is None:
            record = self._record_for_key(station_name, exact, preferred_oid)
            if record is not None:
                point = record["point"]
        if point is None:
            point = group_data["point"]
        if point is None:
            point = UNKNOWN_POINT

        ih = self._linear_value(setup.attrib.get("instrumentHeight"))
        v0_source = setup.attrib.get("orientationAzimuth")
        if v0_source is None:
            v0_source = setup.attrib.get("circleAzimuth")
        v0 = self._direction_to_gon(v0_source)

        station_record = self._record_for_key(station_name, exact, preferred_oid)
        code = station_record["code"] if station_record is not None else None
        if code is None:
            code = group_data["code"]
        if code is None:
            code = self._terrain_code(setup)
        feature = Feature(
            point,
            desc="ST",
            id=pid,
            point_name=station_name,
            st_name=station_name,
            setup_id=setup_id,
            angle_unit="gon",
            dist_unit="meter",
            ih=self._format_number(ih),
            v0=self._format_number(v0),
            hz0=self._format_number(v0),
            code=code,
        )
        return feature

    def _normalise_face(self, angle, z_angle, direct_face):
        if direct_face is not False:
            return angle, z_angle
        if angle is not None:
            angle = (angle + 200.0) % 400.0
        if z_angle is not None:
            z_angle = (400.0 - z_angle) % 400.0
        return angle, z_angle

    @staticmethod
    def _direct_face(raw):
        value = raw.attrib.get("directFace")
        if value is None:
            return None
        return value.strip().lower() in {"true", "1", "yes"}

    @staticmethod
    def _first_attribute(element, names):
        for name in names:
            value = element.attrib.get(name)
            if value is not None:
                return value
        return None

    def _raw_observation_feature(
        self,
        raw,
        setup_id,
        station,
        exact,
        preferred_oid,
        backsights,
        pid,
    ):
        target = raw.find(self._q("TargetPoint"))
        target_name = None
        if target is not None:
            target_name = target.attrib.get("name") or target.attrib.get("desc")
            if not target_name:
                target_name = target.attrib.get("pntRef")
        if not target_name:
            target_name = "point_{}".format(pid)

        purpose = (raw.attrib.get("purpose") or "normal").strip().lower()
        is_reference = purpose in self._REFERENCE_PURPOSES
        matching_backsight = None
        if is_reference:
            matching_backsight = self._matching_backsight(setup_id, target, backsights)

        record = self._resolve_record(
            target, exact, preferred_oid, prefer_control=is_reference
        )
        if matching_backsight is not None and matching_backsight["record"] is not None:
            record = matching_backsight["record"]

        if is_reference and matching_backsight is not None:
            point = matching_backsight["point"]
        elif is_reference and record is not None:
            point = record["point"]
        else:
            point = self._point_from_text(target.text if target is not None else None)
        if point is None and record is not None:
            point = record["point"]
        if point is None:
            point = UNKNOWN_POINT

        raw_angle = raw.attrib.get("horizAngle")
        if raw_angle is not None:
            angle = self._angle_to_gon(raw_angle)
        else:
            angle = self._direction_to_gon(raw.attrib.get("azimuth"))
        z_angle = self._angle_to_gon(raw.attrib.get("zenithAngle"))
        angle, z_angle = self._normalise_face(angle, z_angle, self._direct_face(raw))

        dist = raw.attrib.get("slopeDistance")
        dist_type = "s"
        if dist is None:
            dist = raw.attrib.get("horizDistance")
            dist_type = "h"
        if dist is None:
            dist_type = None
        dist = self._linear_value(dist)
        target_height = self._linear_value(raw.attrib.get("targetHeight"))
        z_distance = self._linear_value(
            self._first_attribute(
                raw,
                ("vertDistance", "elevDifference", "heightDifference"),
            )
        )

        code = record["code"] if record is not None else None
        if code is None:
            code = self._terrain_code(target)
        if code is None:
            code = self._terrain_code(raw)
        desc = "R" if is_reference else "PO"
        kwargs = {
            "point_name": target_name,
            "angle_unit": "gon",
            "z_angle_type": "z" if z_angle is not None else None,
            "dist_unit": "meter",
            "dist_type": dist_type,
            "angle": self._format_number(angle),
            "z_angle": self._format_number(z_angle),
            "dist": self._format_number(dist),
            "z_distance": self._format_number(z_distance),
            "th": self._format_number(target_height),
            "ih": station.properties.get("ih"),
            "st_name": station.properties.get("st_name"),
            "setup_id": setup_id,
            "purpose": purpose,
            "code": code,
        }
        if is_reference:
            kwargs["ref_name"] = target_name
        return Feature(point, desc=desc, id=pid, **kwargs)

    def _reduced_observation_feature(
        self,
        observation,
        setup_id,
        station,
        stations,
        exact,
        preferred_oid,
        backsights,
        pid,
    ):
        target_setup_id = observation.attrib.get("targetSetupID")
        target_station = stations.get(target_setup_id)
        target = observation.find(self._q("TargetPoint"))

        target_name = None
        point = None
        record = None
        if target_station is not None:
            target_name = target_station.properties.get("st_name")
            point = target_station.geometry
        if target is not None:
            target_name = (
                target.attrib.get("name")
                or target.attrib.get("desc")
                or target.attrib.get("pntRef")
                or target_name
            )
            direct_point = self._point_from_text(target.text)
            record = self._resolve_record(target, exact, preferred_oid)
            if direct_point is not None:
                point = direct_point
            elif record is not None and record["point"] is not None:
                point = record["point"]

        point_ref = self._first_attribute(
            observation,
            ("targetPntRef", "targetPointRef", "pntRef"),
        )
        if target_name is None and point_ref:
            target_name = point_ref
        if target_name is None:
            target_name = target_setup_id or observation.attrib.get("name")
        if target_name is None:
            target_name = "point_{}".format(pid)

        if record is None:
            record = self._record_for_key(target_name, exact, preferred_oid)
        if record is None and point_ref:
            record = self._record_for_key(point_ref, exact, preferred_oid)
        if point is None and record is not None:
            point = record["point"]

        purpose = (observation.attrib.get("purpose") or "normal").strip().lower()
        is_reference = purpose in self._REFERENCE_PURPOSES
        matching_backsight = None
        if is_reference:
            for entry in backsights.get(setup_id, []):
                if target_name == entry["name"] or target_name in entry["keys"]:
                    entry["used"] = True
                    matching_backsight = entry
                    break
        if matching_backsight is not None:
            point = matching_backsight["point"] or point
            if matching_backsight["record"] is not None:
                record = matching_backsight["record"]

        if point is None:
            point = UNKNOWN_POINT

        raw_angle = observation.attrib.get("azimuth")
        if raw_angle is not None:
            angle = self._direction_to_gon(raw_angle)
        else:
            angle = self._angle_to_gon(observation.attrib.get("horizAngle"))
        z_angle = self._angle_to_gon(observation.attrib.get("zenithAngle"))

        dist = observation.attrib.get("slopeDistance")
        dist_type = "s"
        if dist is None:
            dist = observation.attrib.get("horizDistance")
            dist_type = "h"
        if dist is None:
            dist_type = None
        dist = self._linear_value(dist)
        z_distance = self._linear_value(
            self._first_attribute(
                observation,
                ("vertDistance", "elevDifference", "heightDifference"),
            )
        )
        target_height = self._linear_value(observation.attrib.get("targetHeight"))

        code = record["code"] if record is not None else None
        if code is None:
            code = self._terrain_code(target)
        if code is None:
            code = self._terrain_code(observation)

        desc = "R" if is_reference else "PO"
        kwargs = {
            "point_name": target_name,
            "angle_unit": "gon",
            "z_angle_type": "z" if z_angle is not None else None,
            "dist_unit": "meter",
            "dist_type": dist_type,
            "angle": self._format_number(angle),
            "z_angle": self._format_number(z_angle),
            "dist": self._format_number(dist),
            "z_distance": self._format_number(z_distance),
            "th": self._format_number(target_height),
            "ih": station.properties.get("ih"),
            "st_name": station.properties.get("st_name"),
            "setup_id": setup_id,
            "target_setup_id": target_setup_id,
            "purpose": purpose,
            "code": code,
        }
        if is_reference:
            kwargs["ref_name"] = target_name
        return Feature(point, desc=desc, id=pid, **kwargs)

    def _observation_features(
        self,
        observations,
        setup_id,
        station,
        stations,
        exact,
        preferred_oid,
        backsights,
        pid,
    ):
        features = []
        for observation in observations:
            tag = observation.tag.rsplit("}", 1)[-1]
            observation_setup_id = observation.attrib.get("setupID") or setup_id
            if tag == "RawObservation":
                feature = self._raw_observation_feature(
                    observation,
                    observation_setup_id,
                    station,
                    exact,
                    preferred_oid,
                    backsights,
                    pid,
                )
            elif tag == "ReducedObservation":
                feature = self._reduced_observation_feature(
                    observation,
                    observation_setup_id,
                    station,
                    stations,
                    exact,
                    preferred_oid,
                    backsights,
                    pid,
                )
            else:
                continue
            features.append(feature)
            pid += 1
        return features, pid

    def _synthetic_backsight_feature(self, entry, setup_id, station, pid):
        backsight = entry["element"]
        name = entry["name"] or "reference_{}".format(pid)
        angle_source = backsight.attrib.get("circle")
        if angle_source is not None:
            angle = self._angle_to_gon(angle_source)
        else:
            angle = self._direction_to_gon(backsight.attrib.get("azimuth"))
        target_height = self._linear_value(backsight.attrib.get("targetHeight"))
        record = entry["record"]
        code = record["code"] if record is not None else None
        if code is None:
            code = self._terrain_code(entry["target"])
        if code is None:
            code = self._terrain_code(backsight)
        return Feature(
            entry["point"] or UNKNOWN_POINT,
            desc="R",
            id=pid,
            point_name=name,
            ref_name=name,
            angle_unit="gon",
            dist_unit="meter",
            angle=self._format_number(angle),
            z_angle=None,
            dist=None,
            dist_type=None,
            th=self._format_number(target_height),
            ih=station.properties.get("ih"),
            st_name=station.properties.get("st_name"),
            setup_id=setup_id,
            purpose="backsight",
            code=code,
        )

    def _standalone_cgpoint_features(
        self,
        records,
        exact,
        survey,
        station_names,
        referenced_target_names,
        pid,
    ):
        features = []
        survey_exists = survey is not None
        for record in records:
            if record["point"] is None or not record["name"]:
                continue
            if exact.get(record["name"]) is not record:
                continue
            if record["name"] in station_names:
                continue

            role = record["role"]
            if survey_exists:
                if role == "measured":
                    continue
                if not role and record["name"] in referenced_target_names:
                    continue
                if not role and record["name"].find("@") >= 0:
                    continue

            features.append(
                Feature(
                    record["point"],
                    desc="PT",
                    id=pid,
                    point_name=record["name"],
                    dist_unit="meter",
                    code=record["code"],
                )
            )
            pid += 1
        return features, pid

    @property
    def raw_line(self):
        """Return LandXML coordinates, stations and observations as features."""

        survey = self.line.find(self._q("Survey"))
        records, exact, preferred_oid = self._build_cgpoint_index(survey)

        setups = [] if survey is None else survey.findall(self._q("InstrumentSetup"))
        station_names = {
            setup.attrib.get("stationName")
            for setup in setups
            if setup.attrib.get("stationName")
        }
        referenced_target_names = set()
        if survey is not None:
            for target in survey.iter(self._q("TargetPoint")):
                referenced_target_names.update(self._point_key_candidates(target))

        points = []
        pid = 0
        standalone, pid = self._standalone_cgpoint_features(
            records,
            exact,
            survey,
            station_names,
            referenced_target_names,
            pid,
        )
        points.extend(standalone)

        if survey is None:
            return points

        groups_by_setup = {}
        observations_by_setup = {}
        for group in survey.findall(self._q("ObservationGroup")):
            default_setup_id = self._setup_id_from_group(group)
            group_setup_ids = []
            if default_setup_id:
                group_setup_ids.append(default_setup_id)
            for child in group:
                tag = child.tag.rsplit("}", 1)[-1]
                if tag not in {"Backsight", "RawObservation", "ReducedObservation"}:
                    continue
                setup_id = child.attrib.get("setupID") or default_setup_id
                if setup_id and setup_id not in group_setup_ids:
                    group_setup_ids.append(setup_id)
                if tag in {"RawObservation", "ReducedObservation"}:
                    observations_by_setup.setdefault(setup_id, []).append(child)
            for setup_id in group_setup_ids:
                groups_by_setup.setdefault(setup_id, []).append(group)

        stations = {}
        station_order = []
        for setup in setups:
            setup_id = setup.attrib.get("id", "")
            station = self._station_feature(
                setup,
                groups_by_setup.get(setup_id, []),
                exact,
                preferred_oid,
                pid,
            )
            pid += 1
            setup_id = station.properties.get("setup_id")
            stations[setup_id] = station
            station_order.append(setup_id)

        backsights = self._backsights(survey, exact, preferred_oid)

        emitted_setups = set()
        for setup_id in station_order:
            station = stations[setup_id]
            points.append(station)
            emitted_setups.add(setup_id)

            observations, pid = self._observation_features(
                observations_by_setup.get(setup_id, []),
                setup_id,
                station,
                stations,
                exact,
                preferred_oid,
                backsights,
                pid,
            )
            points.extend(observations)

            for entry in backsights.get(setup_id, []):
                if not entry["used"]:
                    points.append(
                        self._synthetic_backsight_feature(entry, setup_id, station, pid)
                    )
                    pid += 1

        # Keep observations whose setup is declared only on ObservationGroup or
        # RawObservation, even if the corresponding InstrumentSetup is absent.
        remaining_setup_ids = list(groups_by_setup)
        for setup_id in observations_by_setup:
            if setup_id not in remaining_setup_ids:
                remaining_setup_ids.append(setup_id)
        for setup_id in remaining_setup_ids:
            if setup_id in emitted_setups:
                continue
            groups = groups_by_setup.get(setup_id, [])
            group_data = self._group_station_data(groups, exact, preferred_oid)
            station_name = group_data["name"] or setup_id or "station_{}".format(pid)
            station = Feature(
                group_data["point"] or UNKNOWN_POINT,
                desc="ST",
                id=pid,
                point_name=station_name,
                st_name=station_name,
                setup_id=setup_id,
                angle_unit="gon",
                dist_unit="meter",
                ih=None,
                v0=None,
                hz0=None,
                code=group_data["code"],
            )
            pid += 1
            points.append(station)
            stations[setup_id] = station
            observations, pid = self._observation_features(
                observations_by_setup.get(setup_id, []),
                setup_id,
                station,
                stations,
                exact,
                preferred_oid,
                backsights,
                pid,
            )
            points.extend(observations)

        return points

    @property
    def points(self):
        """Return coordinate features, computing polar points when possible."""

        result = []
        stations = {}
        for feature in self.raw_line:
            if feature.desc in {"PT", "P"}:
                result.append(feature)
                continue
            if feature.desc == "ST":
                stations[feature.properties.get("setup_id")] = feature
                continue
            if feature.desc not in {"PO", "R"}:
                continue

            geometry = feature.geometry
            geometry_z = getattr(geometry, "z", None)
            unknown_z = getattr(UNKNOWN_POINT, "z", None)
            if geometry is not None and not (
                geometry.x == UNKNOWN_POINT.x
                and geometry.y == UNKNOWN_POINT.y
                and geometry_z == unknown_z
            ):
                result.append(
                    Feature(
                        geometry,
                        desc="PT",
                        id=feature.id,
                        point_name=feature.point_name,
                    )
                )
                continue

            station = stations.get(feature.properties.get("setup_id"))
            if station is None:
                continue
            sp = station.geometry
            required = (
                sp is not None,
                feature.properties.get("angle") is not None,
                feature.properties.get("z_angle") is not None,
                feature.properties.get("dist") is not None,
            )
            if not all(required):
                continue
            sp_z = getattr(sp, "z", None)
            if sp_z is None:
                continue

            base_point = BasePoint(
                x=sp.x,
                y=sp.y,
                z=sp_z,
                ih=station.properties.get("ih") or 0.0,
                b_zero_st=station.properties.get("v0") or 0.0,
            )
            polar_point = PolarPoint(
                angle_unit="gon",
                z_angle_type=feature.properties.get("z_angle_type") or "z",
                dist_type=feature.properties.get("dist_type") or "s",
                dist=feature.properties["dist"],
                angle=feature.properties["angle"],
                z_angle=feature.properties["z_angle"],
                th=feature.properties.get("th") or 0.0,
                base_point=base_point,
                pid=feature.id,
                text=feature.point_name,
                coordorder="ENZ",
            )
            result.append(
                Feature(
                    polar_point.to_point(),
                    desc="PT",
                    id=feature.id,
                    point_name=feature.point_name,
                )
            )
        return result
