enok/src/klima/geo.rs

//! Posisjon fra BIM-modellen, og valg av klimasted ut fra den.
//!
//! En IFC-modell kan oppgi hvor bygningen står på to måter: som breddegrad og lengdegrad på
//! `IfcSite`, eller som prosjiserte koordinater gjennom `IfcMapConversion` og
//! `IfcProjectedCRS`. I Norge er det siste normalt EUREF89 UTM sone 32, 33 eller 35. Denne
//! modulen gjør begge om til breddegrad og lengdegrad, slik at nærmeste klimasted kan velges
//! automatisk.

use std::f64::consts::FRAC_PI_2;

/// Store halvakse for GRS80, som EUREF89 bygger på, m.
const A: f64 = 6_378_137.0;
/// Flattrykning for GRS80.
const F: f64 = 1.0 / 298.257_222_101;
/// Målestokkfaktor langs sentralmeridianen i UTM.
const K0: f64 = 0.9996;

#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Posisjon {
    /// Breddegrad, grader nord.
    pub bredde: f64,
    /// Lengdegrad, grader øst.
    pub lengde: f64,
}

impl Posisjon {
    /// Avstand langs jordoverflaten, km.
    pub fn avstand_km(&self, annen: &Posisjon) -> f64 {
        let (b1, b2) = (self.bredde.to_radians(), annen.bredde.to_radians());
        let db = b2 - b1;
        let dl = (annen.lengde - self.lengde).to_radians();
        let h = (db / 2.0).sin().powi(2) + b1.cos() * b2.cos() * (dl / 2.0).sin().powi(2);
        2.0 * 6371.0 * h.sqrt().asin()
    }

    /// Sant for plasseringen Vectorworks setter når prosjektstedet ikke er angitt, og for
    /// nullpunktet. Begge betyder at modellen i praksis ikke oppgir noen posisjon.
    pub fn er_standardplassering(&self) -> bool {
        let washington =
            (self.bredde - 38.8895).abs() < 0.05 && (self.lengde + 77.0353).abs() < 0.05;
        let null = self.bredde.abs() < 1e-6 && self.lengde.abs() < 1e-6;
        washington || null
    }

    /// Grovt: innenfor fastlands-Norge og Svalbard.
    pub fn i_norge(&self) -> bool {
        (57.5..=81.5).contains(&self.bredde) && (4.0..=32.0).contains(&self.lengde)
    }
}

/// Fra Web Mercator (EPSG:3857), som Vectorworks bruker når ingen annen projeksjon er valgt.
pub fn fra_web_mercator(x: f64, y: f64) -> Posisjon {
    Posisjon {
        bredde: (2.0 * (y / A).exp().atan() - FRAC_PI_2).to_degrees(),
        lengde: (x / A).to_degrees(),
    }
}

/// Fra UTM på nordlige halvkule, med GRS80-ellipsoiden som EUREF89 bruker.
///
/// Inversformlene er de vanlige seriene etter Snyder, «Map Projections – A Working Manual»,
/// som gir centimeternøyaktighet innenfor en sone.
pub fn fra_utm(sone: u32, x: f64, y: f64) -> Posisjon {
    let e2 = F * (2.0 - F);
    let ep2 = e2 / (1.0 - e2);
    let e1 = (1.0 - (1.0 - e2).sqrt()) / (1.0 + (1.0 - e2).sqrt());

    let m = y / K0;
    let mu = m / (A * (1.0 - e2 / 4.0 - 3.0 * e2 * e2 / 64.0 - 5.0 * e2.powi(3) / 256.0));
    let fi1 = mu
        + (3.0 * e1 / 2.0 - 27.0 * e1.powi(3) / 32.0) * (2.0 * mu).sin()
        + (21.0 * e1 * e1 / 16.0 - 55.0 * e1.powi(4) / 32.0) * (4.0 * mu).sin()
        + (151.0 * e1.powi(3) / 96.0) * (6.0 * mu).sin()
        + (1097.0 * e1.powi(4) / 512.0) * (8.0 * mu).sin();

    let (s, c, t) = (fi1.sin(), fi1.cos(), fi1.tan());
    let n1 = A / (1.0 - e2 * s * s).sqrt();
    let t1 = t * t;
    let c1 = ep2 * c * c;
    let r1 = A * (1.0 - e2) / (1.0 - e2 * s * s).powf(1.5);
    let d = (x - 500_000.0) / (n1 * K0);

    let bredde = fi1
        - (n1 * t / r1)
            * (d * d / 2.0
                - (5.0 + 3.0 * t1 + 10.0 * c1 - 4.0 * c1 * c1 - 9.0 * ep2) * d.powi(4) / 24.0
                + (61.0 + 90.0 * t1 + 298.0 * c1 + 45.0 * t1 * t1 - 252.0 * ep2 - 3.0 * c1 * c1)
                    * d.powi(6)
                    / 720.0);
    let sentralmeridian = (sone as f64 - 1.0) * 6.0 - 180.0 + 3.0;
    let lengde = sentralmeridian
        + ((d - (1.0 + 2.0 * t1 + c1) * d.powi(3) / 6.0
            + (5.0 - 2.0 * c1 + 28.0 * t1 - 3.0 * c1 * c1 + 8.0 * ep2 + 24.0 * t1 * t1)
                * d.powi(5)
                / 120.0)
            / c)
            .to_degrees();
    Posisjon {
        bredde: bredde.to_degrees(),
        lengde,
    }
}

/// Gjør koordinater i et navngitt koordinatsystem om til posisjon.
///
/// `navn` er slik det står i `IfcProjectedCRS`, for eksempel «EPSG:25832». Koordinatsystemer
/// enok ikke kjenner gir `None`, og da brukes ikke posisjonen.
pub fn fra_koordinatsystem(navn: &str, x: f64, y: f64) -> Option<Posisjon> {
    let stor = navn.to_uppercase();
    let kode: Option<u32> = stor.find("EPSG").and_then(|i| {
        let siffer: String = stor[i + 4..]
            .chars()
            .skip_while(|c| !c.is_ascii_digit())
            .take_while(|c| c.is_ascii_digit())
            .collect();
        siffer.parse().ok()
    });
    match kode {
        Some(3857) | Some(900_913) => Some(fra_web_mercator(x, y)),
        // Geografiske koordinater: østverdien er lengdegrad og nordverdien breddegrad.
        Some(4326) | Some(4258) => Some(Posisjon {
            bredde: y,
            lengde: x,
        }),
        // ETRS89 / EUREF89 UTM sone 28–38.
        Some(k @ 25828..=25838) => Some(fra_utm(k - 25800, x, y)),
        // WGS 84 UTM nord, sone 28–38.
        Some(k @ 32628..=32638) => Some(fra_utm(k - 32600, x, y)),
        // ETRS89 / UTM med NN2000-høyder, som Kartverket bruker: 5972 er sone 32, 5973 sone 33,
        // 5975 sone 35.
        Some(k @ 5972..=5976) => Some(fra_utm(k - 5940, x, y)),
        _ => {
            // Uten EPSG-kode: kjenn igjen «UTM32», «UTM 33N» og lignende.
            let i = stor.find("UTM")?;
            let sone: String = stor[i + 3..]
                .chars()
                .skip_while(|c| !c.is_ascii_digit())
                .take_while(|c| c.is_ascii_digit())
                .collect();
            let sone: u32 = sone.parse().ok()?;
            (1..=60).contains(&sone).then(|| fra_utm(sone, x, y))
        }
    }
}

#[cfg(test)]
mod tester {
    use super::*;

    #[test]
    fn web_mercator_fra_vectorworks_gir_washington() {
        // Verdiene i IfcMapConversion når prosjektstedet ikke er satt i Vectorworks.
        let p = fra_web_mercator(-8_575_524.352_288_59, 4_705_852.257_469_73);
        assert!((p.bredde - 38.8895).abs() < 1e-3, "bredde {}", p.bredde);
        assert!((p.lengde + 77.0352).abs() < 1e-3, "lengde {}", p.lengde);
        assert!(p.er_standardplassering());
        assert!(!p.i_norge());
    }

    /// Referanseverdier fra pyproj (PROJ), transformert fra EPSG:4326.
    #[test]
    fn utm_og_web_mercator_stemmer_med_proj() {
        // (sone, øst, nord, breddegrad, lengdegrad, toleranse i grader)
        //
        // Innenfor sonen og et stykke utenfor er avviket under 1e-5 grader, om lag én meter.
        // Tromsø i sone 32 ligger nesten 10 grader utenfor sonen, der rekkeutviklingen gir
        // om lag 1,3 m avvik. Det er langt innenfor hva valg av klimasted trenger, men viser
        // hvor grensen går.
        let utm = [
            (32, 613_351.291, 6_611_921.980, 59.6300, 11.0100, 1e-5),
            (32, 596_946.981, 6_642_824.393, 59.9115, 10.7336, 1e-5),
            (33, 275_055.615, 6_616_968.071, 59.6300, 11.0100, 1e-5),
            (33, 653_421.188, 7_731_721.083, 69.6492, 18.9553, 1e-5),
            (35, 188_546.727, 7_747_295.774, 69.6492, 18.9553, 1e-5),
            (32, 884_909.622, 7_758_204.184, 69.6492, 18.9553, 1e-4),
        ];
        for (sone, x, y, b, l, toleranse) in utm {
            let p = fra_utm(sone, x, y);
            assert!(
                (p.bredde - b).abs() < toleranse && (p.lengde - l).abs() < toleranse,
                "sone {sone} ({x}, {y}) ga {p:?}, forventet {b}, {l}"
            );
        }
        let p = fra_web_mercator(1_225_627.594, 8_317_818.189);
        assert!(
            (p.bredde - 59.63).abs() < 1e-6 && (p.lengde - 11.01).abs() < 1e-6,
            "{p:?}"
        );
    }

    #[test]
    fn sentralmeridianen_gir_eksakt_lengdegrad() {
        // Østverdi 500 000 ligger på sonens sentralmeridian, 9° øst for sone 32.
        let p = fra_utm(32, 500_000.0, 6_600_000.0);
        assert!((p.lengde - 9.0).abs() < 1e-9, "lengde {}", p.lengde);
        let q = fra_utm(33, 500_000.0, 6_600_000.0);
        assert!((q.lengde - 15.0).abs() < 1e-9, "lengde {}", q.lengde);
    }

    #[test]
    fn ekvator_paa_sentralmeridianen() {
        let p = fra_utm(32, 500_000.0, 0.0);
        assert!(p.bredde.abs() < 1e-9 && (p.lengde - 9.0).abs() < 1e-9);
    }

    #[test]
    fn ukjent_koordinatsystem_gir_ingen_posisjon() {
        assert!(fra_koordinatsystem("EPSG:2154", 650_000.0, 6_860_000.0).is_none());
        assert!(fra_koordinatsystem("lokalt", 1.0, 2.0).is_none());
    }

    #[test]
    fn epsg_koder_tolkes() {
        let a = fra_koordinatsystem("EPSG:25832", 500_000.0, 6_600_000.0).unwrap();
        let b = fra_koordinatsystem("ETRS89 / UTM zone 32N", 500_000.0, 6_600_000.0).unwrap();
        assert!((a.lengde - 9.0).abs() < 1e-9 && (b.lengde - 9.0).abs() < 1e-9);
        let c = fra_koordinatsystem("EPSG:5973", 500_000.0, 6_600_000.0).unwrap();
        assert!(
            (c.lengde - 15.0).abs() < 1e-9,
            "5973 er sone 33, fikk {}",
            c.lengde
        );
    }
}