Tracing Unit Errors from Feet to Metres

import json
import math
from pathlib import Path

import numpy as np

FOOT_INTERNATIONAL = 0.3048
FOOT_US_SURVEY = 1200.0 / 3937.0        # 0.30480060960121924

RATIOS = {
    "feet_read_as_metres": {
        "ratio": 1.0 / FOOT_INTERNATIONAL,
        "value": 3.2808,
        "symptom": "the model is 3.28× too large, or a metric model is 3.28× too "
                   "small",
        "where": "a length in feet passed to code expecting metres",
    },
    "metres_read_as_feet": {
        "ratio": FOOT_INTERNATIONAL,
        "value": 0.3048,
        "symptom": "the model is 3.28× too small",
        "where": "a metric length multiplied by 0.3048 'to convert'",
    },
    "millimetres_read_as_metres": {
        "ratio": 1000.0,
        "value": 1000.0,
        "symptom": "the model is 1000× too large; a 12 m building spans 12 km",
        "where": "an IFC or CAD file in millimetres with no unit handling",
    },
    "centimetres_read_as_metres": {
        "ratio": 100.0,
        "value": 100.0,
        "symptom": "100× too large",
        "where": "a scanner export in centimetres",
    },
    "inches_read_as_metres": {
        "ratio": 1.0 / 0.0254,
        "value": 39.3701,
        "symptom": "39.4× too large",
        "where": "an imperial CAD drawing",
    },
    "square_feet_read_as_square_metres": {
        "ratio": 1.0 / (FOOT_INTERNATIONAL ** 2),
        "value": 10.7639,
        "symptom": "an area is 10.76× out",
        "where": "a declared footprint area against a computed one",
    },
    "cubic_feet_read_as_cubic_metres": {
        "ratio": 1.0 / (FOOT_INTERNATIONAL ** 3),
        "value": 35.3147,
        "symptom": "a volume is 35.3× out",
        "where": "an earthworks or stockpile quantity",
    },
    "us_survey_vs_international_foot": {
        "ratio": FOOT_US_SURVEY / FOOT_INTERNATIONAL,
        "value": 1.000002,
        "symptom": "2 ppm — 1.2 m over a 600 km state plane easting",
        "where": "a US state plane CRS where the foot definition matters",
    },
}

def identify_ratio(observed, expected, tolerance=0.01):
    """Given a measured and an expected value, name the unit error."""
    if expected == 0:
        return {"identified": False, "reason": "expected value is zero"}
    r = observed / expected
    candidates = []
    for name, spec in RATIOS.items():
        for candidate in (spec["value"], 1.0 / spec["value"]):
            if abs(r - candidate) / candidate <= tolerance:
                candidates.append({
                    "error": name,
                    "ratio_observed": round(r, 6),
                    "ratio_expected": round(candidate, 6),
                    "inverted": candidate != spec["value"],
                    "symptom": spec["symptom"],
                    "where": spec["where"],
                })
    return {
        "observed": observed, "expected": expected,
        "ratio": round(r, 6),
        "identified": bool(candidates),
        "candidates": candidates,
        "note": ("no known unit ratio matches; the error is not a unit error"
                 if not candidates else None),
    }

for name, spec in RATIOS.items():
    print(f"{spec['value']:>10.4f}  {name}")

def ifc_units(ifc_path):
    """IFC declares its units; the file is authoritative and often ignored."""
    import ifcopenshell
    import ifcopenshell.util.unit as unit_util

    model = ifcopenshell.open(str(ifc_path))
    assignments = model.by_type("IfcUnitAssignment")
    declared = []
    for assignment in assignments:
        for u in assignment.Units or []:
            if u.is_a("IfcSIUnit"):
                declared.append({
                    "kind": u.UnitType,
                    "type": "SI",
                    "name": u.Name,
                    "prefix": u.Prefix,
                })
            elif u.is_a("IfcConversionBasedUnit"):
                factor = u.ConversionFactor
                declared.append({
                    "kind": u.UnitType,
                    "type": "conversion",
                    "name": u.Name,
                    "factor_to_si": float(factor.ValueComponent.wrappedValue)
                    if factor else None,
                })
    scale = unit_util.calculate_unit_scale(model)
    return {
        "file": Path(ifc_path).name,
        "declared_units": declared,
        "length_scale_to_metres": round(float(scale), 12),
        "is_metres": abs(scale - 1.0) < 1e-9,
        "is_millimetres": abs(scale - 0.001) < 1e-12,
        "is_feet": abs(scale - FOOT_INTERNATIONAL) < 1e-6,
        "advice": ("multiply every coordinate by length_scale_to_metres; never "
                   "assume millimetres because most files are"),
    }

def citygml_units(gml_path):
    """CityGML inherits units from its CRS; there is no separate unit element."""
    import re
    from pyproj import CRS

    text = Path(gml_path).read_text(errors="replace")[:200_000]
    srs = re.search(r'srsName="([^"]+)"', text)
    if not srs:
        return {"declared": False,
                "advice": "no srsName found; the unit is unknown and must be "
                          "supplied out of band"}
    name = srs.group(1)
    try:
        crs = CRS.from_user_input(name.split("::")[-1] if "::" in name else name)
    except Exception as exc:
        return {"declared": True, "srsName": name, "parsed": False,
                "error": repr(exc)[:120]}
    axis = crs.axis_info[0] if crs.axis_info else None
    return {
        "declared": True,
        "srsName": name,
        "epsg": crs.to_epsg(),
        "crs_name": crs.name,
        "unit": axis.unit_name if axis else None,
        "unit_conversion_factor": axis.unit_conversion_factor if axis else None,
        "is_metre": bool(axis and axis.unit_name in ("metre", "meter")),
        "advice": ("CityGML has no unit element — the CRS decides, so a file in a "
                   "foot-based state plane CRS has foot coordinates"),
    }

def gltf_units(glb_path):
    """glTF is defined to be in metres. There is no unit field, and that is the trap."""
    import struct
    data = Path(glb_path).read_bytes()
    json_len = struct.unpack_from("<I", data, 12)[0]
    doc = json.loads(data[20:20 + json_len].decode("utf-8"))
    extras = doc.get("asset", {}).get("extras", {})
    return {
        "file": Path(glb_path).name,
        "generator": doc.get("asset", {}).get("generator"),
        "declared_unit": "metres by specification",
        "extras_unit_hint": extras.get("unit") or extras.get("units"),
        "node_scales": [n.get("scale") for n in doc.get("nodes", [])
                        if n.get("scale")][:4],
        "advice": ("glTF has no unit field: the specification says metres. A node "
                   "scale of 0.3048 or 0.001 is a converter compensating for a "
                   "non-metric source, and it should have been baked in"),
    }

KNOWN_DIMENSIONS = {
    "door_height_m": (1.9, 2.4),
    "storey_height_m": (2.2, 6.0),
    "residential_building_height_m": (2.5, 40.0),
    "parcel_frontage_m": (4.0, 200.0),
    "road_lane_width_m": (2.5, 4.2),
    "kerb_height_m": (0.08, 0.25),
    "city_extent_km": (2.0, 60.0),
}

def plausibility_report(measurements, known=KNOWN_DIMENSIONS):
    """measurements: {'storey_height_m': 7.87, ...}"""
    rows = []
    for name, value in measurements.items():
        bounds = known.get(name)
        if bounds is None:
            rows.append({"measurement": name, "value": value,
                         "status": "no reference range"})
            continue
        lo, hi = bounds
        if lo <= value <= hi:
            rows.append({"measurement": name, "value": value,
                         "range": bounds, "status": "plausible"})
            continue
        midpoint = (lo + hi) / 2.0
        diagnosis = identify_ratio(value, midpoint, tolerance=0.15)
        rows.append({
            "measurement": name, "value": value, "range": bounds,
            "status": "IMPLAUSIBLE",
            "ratio_to_midpoint": round(value / midpoint, 4),
            "likely_error": (diagnosis["candidates"][0]["error"]
                             if diagnosis["candidates"] else "not a known unit ratio"),
        })
    bad = [r for r in rows if r["status"] == "IMPLAUSIBLE"]
    ratios = [r["ratio_to_midpoint"] for r in bad if "ratio_to_midpoint" in r]
    consistent = (len(ratios) > 1
                  and max(ratios) / max(min(ratios), 1e-9) < 1.3)
    return {
        "rows": rows,
        "implausible": len(bad),
        "ratios": ratios,
        "one_consistent_ratio": consistent,
        "verdict": (f"every implausible measurement is out by about "
                    f"{np.median(ratios):.2f}× — a single unit error upstream"
                    if consistent else
                    "measurements are out by different factors — not one unit error"
                    if bad else "all measurements plausible"),
    }

def measure_from_mesh(mesh_path):
    import trimesh
    mesh = trimesh.load(mesh_path, process=False, force="mesh")
    extents = mesh.extents
    return {
        "residential_building_height_m": round(float(extents[2]), 3),
        "parcel_frontage_m": round(float(max(extents[0], extents[1])), 3),
        "_bbox_m": [round(float(v), 3) for v in extents],
        "_volume_m3": round(float(mesh.volume), 3) if mesh.is_watertight else None,
        "_area_m2": round(float(mesh.area), 3),
    }

def stage_measurement(path, label, kind="mesh"):
    if kind == "mesh":
        m = measure_from_mesh(path)
    elif kind == "ifc":
        u = ifc_units(path)
        m = {"_length_scale": u["length_scale_to_metres"],
             "_declared": u["declared_units"][:2]}
    elif kind == "gml":
        u = citygml_units(path)
        m = {"_srs": u.get("srsName"), "_unit": u.get("unit"),
             "_factor": u.get("unit_conversion_factor")}
    elif kind == "glb":
        u = gltf_units(path)
        m = {"_generator": u["generator"], "_node_scales": u["node_scales"],
             "_hint": u["extras_unit_hint"]}
    else:
        m = {}
    return {"stage": label, "path": str(path), "kind": kind, **m}

def bisect_units(stage_specs, reference_measurement, reference_value,
                 tolerance=0.02):
    """stage_specs: [(label, path, kind), ...] in pipeline order."""
    probes = []
    for label, path, kind in stage_specs:
        probe = stage_measurement(path, label, kind)
        measured = probe.get(reference_measurement)
        if measured is not None:
            ratio = measured / max(reference_value, 1e-12)
            probe["ratio_to_reference"] = round(ratio, 6)
            probe["correct"] = abs(ratio - 1.0) <= tolerance
            if not probe["correct"]:
                diag = identify_ratio(measured, reference_value)
                probe["likely_error"] = (diag["candidates"][0]["error"]
                                         if diag["candidates"] else "unknown ratio")
        probes.append(probe)

    measured_probes = [p for p in probes if "correct" in p]
    first_bad = next((i for i, p in enumerate(measured_probes)
                      if not p["correct"]), None)
    if first_bad is None:
        return {"probes": probes, "culprit": None,
                "verdict": "every measurable stage is correct; the error is in a "
                           "stage this reference cannot measure"}
    culprit = measured_probes[first_bad]
    previous = measured_probes[first_bad - 1] if first_bad > 0 else None
    return {
        "probes": probes,
        "culprit_stage": culprit["stage"],
        "ratio": culprit["ratio_to_reference"],
        "likely_error": culprit.get("likely_error"),
        "last_correct_stage": previous["stage"] if previous else None,
        "verdict": (f"'{culprit['stage']}' introduced a factor of "
                    f"{culprit['ratio_to_reference']} "
                    f"({culprit.get('likely_error')})"),
    }

def cross_field_unit_check(records, geometry_area_key="computed_area_m2",
                           declared_area_key="footprint_area_m2",
                           height_key="measured_height_m",
                           storeys_key="storeys_above_ground"):
    """A unit error often lands on one attribute and not the geometry."""
    rows = []
    for r in records:
        computed = r.get(geometry_area_key)
        declared = r.get(declared_area_key)
        height = r.get(height_key)
        storeys = r.get(storeys_key)

        entry = {"id": r.get("register_id", "?")}
        if computed and declared:
            entry["area_ratio"] = round(declared / max(computed, 1e-9), 4)
        if height and storeys:
            entry["per_storey_m"] = round(height / max(storeys, 1e-9), 3)
        rows.append(entry)

    area_ratios = [r["area_ratio"] for r in rows if "area_ratio" in r]
    per_storey = [r["per_storey_m"] for r in rows if "per_storey_m" in r]

    findings = []
    if area_ratios:
        median = float(np.median(area_ratios))
        diag = identify_ratio(median, 1.0, tolerance=0.03)
        if abs(median - 1.0) > 0.03:
            findings.append({
                "field": declared_area_key,
                "median_ratio": round(median, 4),
                "likely_error": (diag["candidates"][0]["error"]
                                 if diag["candidates"] else "unknown"),
                "note": "the declared area disagrees with the geometry, so only the "
                        "attribute was converted",
            })
    if per_storey:
        median = float(np.median(per_storey))
        if not 2.2 <= median <= 6.0:
            diag = identify_ratio(median, 3.0, tolerance=0.2)
            findings.append({
                "field": height_key,
                "median_per_storey_m": round(median, 3),
                "likely_error": (diag["candidates"][0]["error"]
                                 if diag["candidates"] else "unknown"),
                "note": "height and storey count disagree, so the height alone was "
                        "converted or not converted",
            })
    return {
        "records": len(rows),
        "median_area_ratio": round(float(np.median(area_ratios)), 4)
        if area_ratios else None,
        "median_per_storey_m": round(float(np.median(per_storey)), 3)
        if per_storey else None,
        "findings": findings,
        "clean": not findings,
    }

def convert_to_metres(values, from_unit, us_survey_foot=False):
    factors = {
        "m": 1.0, "metre": 1.0, "meter": 1.0,
        "mm": 0.001, "millimetre": 0.001,
        "cm": 0.01, "centimetre": 0.01,
        "km": 1000.0,
        "ft": FOOT_US_SURVEY if us_survey_foot else FOOT_INTERNATIONAL,
        "foot": FOOT_US_SURVEY if us_survey_foot else FOOT_INTERNATIONAL,
        "in": 0.0254, "inch": 0.0254,
        "yd": 0.9144, "yard": 0.9144,
    }
    key = from_unit.strip().lower()
    if key not in factors:
        raise ValueError(f"unknown unit {from_unit!r}; declare it explicitly")
    return np.asarray(values, dtype=np.float64) * factors[key], {
        "from_unit": key,
        "factor": factors[key],
        "us_survey_foot": us_survey_foot and key in ("ft", "foot"),
    }

class Metres:
    """A boundary type: once inside the pipeline, everything is metres."""
    __slots__ = ("values", "provenance")

    def __init__(self, values, provenance):
        self.values = np.asarray(values, dtype=np.float64)
        self.provenance = provenance

    @classmethod
    def from_unit(cls, values, unit, us_survey_foot=False):
        converted, prov = convert_to_metres(values, unit,
                                            us_survey_foot=us_survey_foot)
        return cls(converted, prov)

    def assert_plausible(self, expected_range, label="length"):
        lo, hi = expected_range
        span = float(np.ptp(self.values)) if self.values.size > 1 else \
            float(self.values.max())
        if not lo <= span <= hi:
            diag = identify_ratio(span, (lo + hi) / 2.0)
            raise AssertionError(
                f"{label} span {span:.3f} m is outside {expected_range}; "
                f"likely {diag['candidates'][0]['error'] if diag['candidates'] else 'not a unit error'}"
                f" (converted from {self.provenance['from_unit']})")
        return {"label": label, "span_m": round(span, 3),
                "range": expected_range, "ok": True,
                "provenance": self.provenance}

def ingest_ifc(ifc_path, expected_height_range=(2.5, 400.0)):
    """Convert once, at the boundary, and assert before anything else runs."""
    import ifcopenshell
    import ifcopenshell.util.unit as unit_util
    import numpy as np

    model = ifcopenshell.open(str(ifc_path))
    scale = float(unit_util.calculate_unit_scale(model))
    if not 0.0001 <= scale <= 1000.0:
        raise ValueError(f"implausible unit scale {scale} in {ifc_path}")

    # Collect a height from the model's own extent for the assertion.
    import ifcopenshell.geom
    settings = ifcopenshell.geom.settings()
    settings.set(settings.USE_WORLD_COORDS, True)
    zs = []
    for product in model.by_type("IfcBuildingElement")[:200]:
        if not product.Representation:
            continue
        try:
            shape = ifcopenshell.geom.create_shape(settings, product)
        except Exception:
            continue
        verts = np.asarray(shape.geometry.verts, dtype=np.float64).reshape(-1, 3)
        if verts.size:
            zs.extend(verts[:, 2].tolist())
    if not zs:
        raise ValueError("no geometry read; cannot assert the unit")

    heights = Metres.from_unit(np.asarray(zs), "m")   # geometry is already scaled
    heights.values *= 1.0                              # explicit: no double scaling
    check = heights.assert_plausible(expected_height_range, label="building height")
    return {"ifc": Path(ifc_path).name, "declared_scale_to_metres": scale,
            "assertion": check,
            "note": "ifcopenshell's create_shape already applies the unit scale; "
                    "multiplying again is the commonest double-conversion bug"}

    3.2808  feet_read_as_metres
    0.3048  metres_read_as_feet
 1000.0000  millimetres_read_as_metres
  100.0000  centimetres_read_as_metres
   39.3701  inches_read_as_metres
   10.7639  square_feet_read_as_square_metres
   35.3147  cubic_feet_read_as_cubic_metres
    1.0000  us_survey_vs_international_foot
{
  "file": "building_a.ifc",
  "declared_units": [{"kind": "LENGTHUNIT", "type": "SI", "name": "METRE",
                      "prefix": "MILLI"}],
  "length_scale_to_metres": 0.001,
  "is_metres": false, "is_millimetres": true, "is_feet": false
}
{
  "rows": [
    {"measurement": "residential_building_height_m", "value": 46.59,
     "range": [2.5, 40.0], "status": "IMPLAUSIBLE",
     "ratio_to_midpoint": 2.1925, "likely_error": "not a known unit ratio"},
    {"measurement": "parcel_frontage_m", "value": 118.4,
     "range": [4.0, 200.0], "status": "plausible"}
  ],
  "implausible": 1,
  "verdict": "measurements are out by different factors — not one unit error"
}
{
  "observed": 46.59, "expected": 14.2, "ratio": 3.280282,
  "identified": true,
  "candidates": [
    {"error": "feet_read_as_metres", "ratio_observed": 3.280282,
     "ratio_expected": 3.2808, "inverted": false,
     "symptom": "the model is 3.28× too large, or a metric model is 3.28× too small",
     "where": "a length in feet passed to code expecting metres"}
  ]
}
{
  "culprit_stage": "crs_convert",
  "ratio": 3.280282,
  "likely_error": "feet_read_as_metres",
  "last_correct_stage": "simplify",
  "verdict": "'crs_convert' introduced a factor of 3.280282 (feet_read_as_metres)"
}

def double_conversion_check(stage_specs, reference_measurement, reference_value,
                            tolerance=0.02):
    """After fixing, no stage should apply the factor a second time."""
    probes = []
    for label, path, kind in stage_specs:
        probe = stage_measurement(path, label, kind)
        measured = probe.get(reference_measurement)
        if measured is None:
            probes.append({**probe, "measurable": False})
            continue
        ratio = measured / max(reference_value, 1e-12)
        probes.append({**probe, "measurable": True,
                       "ratio": round(ratio, 6),
                       "correct": abs(ratio - 1.0) <= tolerance})

    measurable = [p for p in probes if p.get("measurable")]
    ratios = [p["ratio"] for p in measurable]
    changes = []
    for a, b in zip(measurable, measurable[1:]):
        step = b["ratio"] / max(a["ratio"], 1e-12)
        if abs(step - 1.0) > tolerance:
            diag = identify_ratio(step, 1.0, tolerance=0.05)
            changes.append({
                "from": a["stage"], "to": b["stage"],
                "step_ratio": round(step, 6),
                "likely": (diag["candidates"][0]["error"]
                           if diag["candidates"] else "unknown"),
            })
    return {
        "stages_measured": len(measurable),
        "ratios": [round(r, 4) for r in ratios],
        "all_correct": all(p["correct"] for p in measurable),
        "unit_changes_between_stages": changes,
        "double_conversion": len(changes) > 1,
        "verdict": ("one conversion at one boundary — correct"
                    if len(changes) <= 1 and all(p["correct"] for p in measurable)
                    else f"{len(changes)} unit changes detected; a conversion is "
                         f"being applied more than once"),
    }

def survey_foot_check(easting_ft, northing_ft, epsg_code, tolerance_m=0.05):
    """For a US state plane CRS in feet, which foot definition was intended?"""
    from pyproj import CRS

    crs = CRS.from_epsg(epsg_code)
    axis = crs.axis_info[0] if crs.axis_info else None
    declared = axis.unit_name if axis else "unknown"
    factor = axis.unit_conversion_factor if axis else None

    as_international = np.array([easting_ft, northing_ft]) * FOOT_INTERNATIONAL
    as_survey = np.array([easting_ft, northing_ft]) * FOOT_US_SURVEY
    difference_m = float(np.linalg.norm(as_survey - as_international))

    return {
        "epsg": epsg_code,
        "crs_name": crs.name,
        "declared_unit": declared,
        "declared_factor": factor,
        "matches_international": factor is not None
        and abs(factor - FOOT_INTERNATIONAL) < 1e-9,
        "matches_us_survey": factor is not None
        and abs(factor - FOOT_US_SURVEY) < 1e-12,
        "position_difference_m": round(difference_m, 4),
        "matters": difference_m > tolerance_m,
        "advice": ("use the factor the CRS declares; the US survey foot was "
                   "deprecated at the end of 2022 but remains in existing data, "
                   "and mixing the two is a 2 ppm error"
                   if difference_m > tolerance_m
                   else "the difference is below tolerance at this coordinate "
                        "magnitude"),
    }

print(json.dumps(survey_foot_check(2_184_120.4, 618_402.8, 2229), indent=2))