Source code for openworld_radio_twin.geodesy

from dataclasses import dataclass

from pyproj import Transformer
from pyproj.enums import TransformDirection




@dataclass(frozen=True)
class LocalFrame:
    """WGS84 Earth-centered to local East-North-Up transformation."""

    origin_longitude: float
    origin_latitude: float
    transformer: Transformer



    @classmethod
    def at(cls, longitude: float, latitude: float) -> "LocalFrame":
        pipeline = (
            "+proj=pipeline "
            "+step +proj=cart +ellps=WGS84 "
            f"+step +proj=topocentric +ellps=WGS84 +lon_0={longitude:.12f} "
            f"+lat_0={latitude:.12f} +h_0=0"
        )
        return cls(longitude, latitude, Transformer.from_pipeline(pipeline))




    def to_local(
        self, longitude: float, latitude: float, altitude_m: float = 0
    ) -> tuple[float, float, float]:
        east, north, up = self.transformer.transform(longitude, latitude, altitude_m)
        return float(east), float(north), float(up)




    def to_geographic(
        self, east_m: float, north_m: float, up_m: float = 0
    ) -> tuple[float, float, float]:
        longitude, latitude, altitude = self.transformer.transform(
            east_m,
            north_m,
            up_m,
            direction=TransformDirection.INVERSE,
        )
        return float(longitude), float(latitude), float(altitude)




    def corners(
        self, west_m: float, south_m: float, east_m: float, north_m: float
    ) -> list[list[float]]:
        local_corners = [
            (west_m, north_m),
            (east_m, north_m),
            (east_m, south_m),
            (west_m, south_m),
        ]
        return [list(self.to_geographic(east, north)[:2]) for east, north in local_corners]