Unit Conversion Pipelines for CAD, GIS & BIM Interoperability

# pint>=0.23, ezdxf>=1.2, ifcopenshell
from __future__ import annotations
from typing import Optional
import logging
import ezdxf
import ifcopenshell

log = logging.getLogger(__name__)

# DXF $INSUNITS → unit string mapping
INSUNITS_MAP: dict[int, str] = {
    1: "inch", 2: "foot", 4: "millimeter",
    5: "centimeter", 6: "meter", 7: "kilometer",
}

def extract_dxf_unit(path: str) -> str:
    doc = ezdxf.readfile(path)
    insunits: int = doc.header.get("$INSUNITS", 0)
    if insunits == 0 or insunits not in INSUNITS_MAP:
        # Heuristic: sample bounding box of modelspace
        unit = _heuristic_dxf_unit(doc)
        log.warning("$INSUNITS absent/unitless in %s; inferred %s", path, unit)
        return unit
    return INSUNITS_MAP[insunits]

def _heuristic_dxf_unit(doc: ezdxf.document.Drawing) -> str:
    """Infer unit from coordinate magnitude of modelspace entities."""
    msp = doc.modelspace()
    extents = [e.dxf.insert for e in msp if hasattr(e.dxf, "insert")]
    if not extents:
        return "millimeter"  # safe default for architectural CAD
    max_coord = max(abs(p.x) + abs(p.y) for p in extents)
    if max_coord > 50_000:
        return "millimeter"
    if max_coord > 500:
        return "meter"
    return "foot"

def extract_ifc_unit(path: str) -> str:
    model = ifcopenshell.open(path)
    project = model.by_type("IfcProject")[0]
    for unit_assignment in project.UnitsInContext.Units:
        if getattr(unit_assignment, "UnitType", None) == "LENGTHUNIT":
            prefix = getattr(unit_assignment, "Prefix", None) or ""
            name = getattr(unit_assignment, "Name", "METRE")
            return f"{prefix.lower()}{name.lower()}".replace("metre", "meter")
    log.warning("No length unit found in IfcUnitAssignment; defaulting to meter")
    return "meter"

# pint>=0.23, numpy>=1.26
from __future__ import annotations
import numpy as np
import pint
from dataclasses import dataclass

ureg = pint.UnitRegistry()
Q_ = ureg.Quantity

@dataclass(frozen=True)
class ConversionConfig:
    source_unit: str        # e.g. "millimeter", "foot"
    target_unit: str = "meter"
    tolerance_m: float = 1e-6     # geometric epsilon in target unit
    max_scale_factor: float = 1e4  # sanity guard

class UnitConverter:
    def __init__(self, config: ConversionConfig) -> None:
        self.config = config
        self.scale_factor: float = (
            Q_(1, config.source_unit).to(config.target_unit).magnitude
        )
        self._validate()

    def _validate(self) -> None:
        if not (0 < self.scale_factor <= self.config.max_scale_factor):
            raise ValueError(
                f"Scale factor {self.scale_factor:.6g} out of safe range "
                f"({self.config.source_unit} → {self.config.target_unit})"
            )

    def scale_coordinates(self, coords: np.ndarray) -> np.ndarray:
        """Vectorized Nx2 or Nx3 coordinate scaling. Input dtype: float64."""
        if coords.ndim != 2 or coords.shape[1] not in (2, 3):
            raise ValueError("coords must be Nx2 or Nx3 float64 array")
        return coords * self.scale_factor

    def scale_attribute(self, values: np.ndarray, attr_unit: str) -> np.ndarray:
        """Scale a non-geometric attribute with its own declared unit."""
        factor = Q_(1, attr_unit).to(self.config.target_unit).magnitude
        return values * factor

# numpy>=1.26, ezdxf>=1.2
import numpy as np
from typing import Iterator
import ezdxf
from ezdxf.entities import DXFGraphic

def iter_entity_coords(
    msp: ezdxf.layouts.Modelspace,
) -> Iterator[tuple[str, np.ndarray]]:
    """Yield (entity_handle, Nx3 coordinate array) for each geometry entity."""
    for entity in msp:
        handle = entity.dxf.handle
        if entity.dxftype() == "LINE":
            pts = np.array([
                [entity.dxf.start.x, entity.dxf.start.y, entity.dxf.start.z],
                [entity.dxf.end.x, entity.dxf.end.y, entity.dxf.end.z],
            ], dtype=np.float64)
            yield handle, pts
        elif entity.dxftype() == "LWPOLYLINE":
            pts = np.array(
                [[p[0], p[1], 0.0] for p in entity.get_points()],
                dtype=np.float64,
            )
            yield handle, pts
        # Extend for ARC, CIRCLE, INSERT (block references), SPLINE, etc.

def convert_modelspace(
    doc: ezdxf.document.Drawing,
    converter: UnitConverter,
) -> dict[str, np.ndarray]:
    """Return {handle: scaled_coords} for the entire modelspace."""
    msp = doc.modelspace()
    result: dict[str, np.ndarray] = {}
    for handle, coords in iter_entity_coords(msp):
        result[handle] = converter.scale_coordinates(coords)
    return result

# shapely>=2.0, numpy>=1.26
from __future__ import annotations
import numpy as np
from shapely.geometry import Polygon
from dataclasses import dataclass

@dataclass
class ValidationResult:
    handle: str
    passed: bool
    message: str

def validate_polygon_closure(
    handle: str,
    coords: np.ndarray,
    epsilon: float = 1e-6,
) -> ValidationResult:
    """Check that a polygon's first and last vertex are within epsilon."""
    if coords.shape[0] < 3:
        return ValidationResult(handle, False, "Too few vertices for polygon")
    gap = np.linalg.norm(coords[0] - coords[-1])
    if gap > epsilon:
        return ValidationResult(
            handle, False,
            f"Open polygon: gap={gap:.3e} m (epsilon={epsilon:.3e} m)"
        )
    poly = Polygon(coords[:, :2])
    if not poly.is_valid:
        return ValidationResult(handle, False, f"Invalid topology: {poly.explain_validity()}")
    return ValidationResult(handle, True, "OK")

def validate_bounding_box(
    original_coords: np.ndarray,
    scaled_coords: np.ndarray,
    expected_scale: float,
    tolerance: float = 1e-4,
) -> bool:
    """Assert that the bounding-box diagonal scaled by expected_scale matches."""
    orig_diag = np.linalg.norm(original_coords.max(0) - original_coords.min(0))
    scaled_diag = np.linalg.norm(scaled_coords.max(0) - scaled_coords.min(0))
    if orig_diag == 0:
        return True
    ratio = scaled_diag / orig_diag
    return abs(ratio - expected_scale) < tolerance

# pyproj>=3.6
from pyproj import CRS, Transformer
import numpy as np

def guard_geographic_crs(epsg: int) -> bool:
    """Return True if the CRS is geographic (degree-based), not projected."""
    crs = CRS.from_epsg(epsg)
    return crs.is_geographic

# pytest>=7.0, numpy>=1.26, pint>=0.23
import numpy as np
import pytest
from your_pipeline import UnitConverter, ConversionConfig  # adjust import

def test_millimeter_to_meter_round_trip() -> None:
    cfg = ConversionConfig(source_unit="millimeter", target_unit="meter")
    conv = UnitConverter(cfg)
    original = np.array([[12500.0, 8000.0, 3000.0]], dtype=np.float64)

    scaled = conv.scale_coordinates(original)
    assert scaled.shape == original.shape
    np.testing.assert_allclose(scaled, [[12.5, 8.0, 3.0]], rtol=1e-9)

    # Reverse pass
    inv_cfg = ConversionConfig(source_unit="meter", target_unit="millimeter")
    inv_conv = UnitConverter(inv_cfg)
    reconstructed = inv_conv.scale_coordinates(scaled)
    np.testing.assert_allclose(reconstructed, original, rtol=1e-9)

def test_invalid_scale_factor_raises() -> None:
    cfg = ConversionConfig(
        source_unit="millimeter", target_unit="meter", max_scale_factor=1e-3
    )
    with pytest.raises(ValueError, match="out of safe range"):
        UnitConverter(cfg)

def test_heuristic_unit_detection_large_coords() -> None:
    """Coordinates > 50 000 should infer millimeters."""
    from your_pipeline import _heuristic_from_max_coord
    assert _heuristic_from_max_coord(75_000) == "millimeter"
    assert _heuristic_from_max_coord(800) == "meter"

def test_attribute_unit_independence() -> None:
    cfg = ConversionConfig(source_unit="foot", target_unit="meter")
    conv = UnitConverter(cfg)
    # Geometry: feet → meters
    geo = np.array([[10.0, 0.0, 0.0]], dtype=np.float64)
    scaled_geo = conv.scale_coordinates(geo)
    np.testing.assert_allclose(scaled_geo[0, 0], 3.048, rtol=1e-6)
    # Attribute declared in millimeters, not feet
    attr = np.array([500.0])  # 500 mm pipe diameter
    scaled_attr = conv.scale_attribute(attr, attr_unit="millimeter")
    np.testing.assert_allclose(scaled_attr[0], 0.5, rtol=1e-6)

# Batch extraction pattern (numpy>=1.26, ezdxf>=1.2)
all_coords = np.vstack([coords for _, coords in iter_entity_coords(msp)])
all_scaled = all_coords * converter.scale_factor
# Split back by entity length using np.split and saved index offsets

# geopandas>=0.14, shapely>=2.0
import geopandas as gpd
from shapely.affinity import scale as shp_scale

CHUNK = 50_000

def convert_gdf_chunks(path: str, factor: float):
    for chunk in gpd.read_file(path, rows=slice(0, None, CHUNK)):
        chunk["geometry"] = chunk["geometry"].apply(
            lambda g: shp_scale(g, xfact=factor, yfact=factor, zfact=factor, origin=(0, 0, 0))
        )
        yield chunk