geo_omerc/
lib.rs

1#![doc = include_str!("../README.md")]
2
3use geo::{Coord, Geometry, MapCoords, Point, coord};
4use geodesy::Direction::{Fwd, Inv};
5use geodesy::Error as GeodesyError;
6use geodesy::coord::{Coor2D, CoordinateTuple};
7use geodesy::ctx::{Context, Minimal, OpHandle};
8
9/// Catch-all for all definition or transformation errors.
10#[derive(Debug)]
11pub struct OMercError;
12
13impl From<GeodesyError> for OMercError {
14    fn from(_: GeodesyError) -> Self {
15        OMercError
16    }
17}
18
19/// An utility to convert `geo::Geometry` types back and forward between
20/// geodesic degrees and cartesian meters.
21pub struct OMercTransformer {
22    context: Minimal,
23    operation: OpHandle,
24}
25
26impl OMercTransformer {
27    /// Creates a new transformer with a geodesic `anchor_point` as the new center
28    /// point (0, 0) of a cartesian plane in meters, projected via Oblique Mercator.
29    pub fn new(anchor_point: &Point) -> Result<Self, OMercError> {
30        let mut context = Minimal::new();
31        let operation = context.op(&format!(
32            "omerc ellps=WGS84 latc={} lonc={} k_0=1.0 alpha=1e-06",
33            anchor_point.y(),
34            anchor_point.x()
35        ))?;
36        Ok(Self { context, operation })
37    }
38
39    /// Returns a new `geo::Geometry` with its coordinates transformed
40    /// from geodesic (in degrees) to cartesian (in meters).
41    pub fn to_cartesian(&self, geometry: &Geometry) -> Result<Geometry, OMercError> {
42        Ok(
43            geometry.try_map_coords(|c: Coord<f64>| -> Result<_, GeodesyError> {
44                let mut data = [Coor2D::geo(c.y, c.x)];
45                self.context.apply(self.operation, Fwd, &mut data)?;
46                Ok(coord! {x: data[0].x(), y: data[0].y()})
47            })?,
48        )
49    }
50
51    /// Returns a new `geo::Geometry` with its coordinates transformed
52    /// from  cartesian (in meters) to geodesic (in degrees).
53    pub fn to_geodesic(&self, geometry: &Geometry) -> Result<Geometry, OMercError> {
54        Ok(
55            geometry.try_map_coords(|c: Coord<f64>| -> Result<_, GeodesyError> {
56                let mut data = [Coor2D::raw(c.x, c.y)];
57                self.context.apply(self.operation, Inv, &mut data)?;
58                let xy = data[0].xy_to_degrees();
59                Ok(coord! {x: xy.0, y: xy.1})
60            })?,
61        )
62    }
63}
64
65#[cfg(test)]
66mod tests {
67    use super::*;
68    use assert_approx_eq::assert_approx_eq;
69    use geo::{Area, Centroid, Distance, Euclidean, polygon};
70
71    #[test]
72    fn roundtrip() {
73        let praca_se = Point::new(-46.632913, -23.550617);
74        let masp = Point::new(-46.655945, -23.561442);
75        let transformer = OMercTransformer::new(&praca_se).unwrap();
76        let c_masp = transformer.to_cartesian(&masp.into()).unwrap();
77        let g_masp: Point = transformer
78            .to_geodesic(&c_masp.into())
79            .unwrap()
80            .try_into()
81            .unwrap();
82        assert_approx_eq!(g_masp.x(), masp.x());
83        assert_approx_eq!(g_masp.y(), masp.y());
84    }
85
86    #[test]
87    fn distance_in_meters() {
88        let center_of_berlin = Point::new(13.409408, 52.520842);
89        let brandenburg_gate = Point::new(13.377701, 52.516262);
90        let transformer = OMercTransformer::new(&center_of_berlin).unwrap();
91        let c_center_of_berlin = transformer.to_cartesian(&center_of_berlin.into()).unwrap();
92        let c_brandenburg_gate = transformer.to_cartesian(&brandenburg_gate.into()).unwrap();
93        let distance = Euclidean.distance(&c_center_of_berlin, &c_brandenburg_gate);
94        assert_approx_eq!(distance, 2211.840487); // meters.
95    }
96
97    #[test]
98    fn area_of_polygon() {
99        let museum_island = polygon![
100            (x: 13.393639, y: 52.522081),
101            (x: 13.397824, y: 52.519285),
102            (x: 13.403055, y: 52.512100),
103            (x: 13.405072, y: 52.511486),
104            (x: 13.411649, y: 52.513873),
105            (x: 13.406418, y: 52.514692),
106            (x: 13.401411, y: 52.519580),
107            (x: 13.398160, y: 52.521513),
108            (x: 13.393639, y: 52.522081),
109        ];
110        let transformer = OMercTransformer::new(&museum_island.centroid().unwrap()).unwrap();
111        let c_museum_island = transformer.to_cartesian(&museum_island.into()).unwrap();
112        assert_approx_eq!(c_museum_island.unsigned_area(), 377902.962989); // square meters.
113    }
114}
geo_omerc/lib.rs

geo_omerc/
lib.rs
