diff --git a/js/src/algorithm/geo/geodesic_area.rs b/js/src/algorithm/geo/geodesic_area.rs new file mode 100644 index 000000000..dbad228ab --- /dev/null +++ b/js/src/algorithm/geo/geodesic_area.rs @@ -0,0 +1,142 @@ +use crate::array::*; +use wasm_bindgen::prelude::*; + +macro_rules! impl_geodesic_area { + ($struct_name:ident) => { + #[wasm_bindgen] + impl $struct_name { + /// Determine the area of a geometry on an ellipsoidal model of the earth. + /// + /// This uses the geodesic measurement methods given by [Karney (2013)]. + /// + /// # Assumptions + /// - Polygons are assumed to be wound in a counter-clockwise direction + /// for the exterior ring and a clockwise direction for interior rings. + /// This is the standard winding for geometries that follow the Simple Feature standard. + /// Alternative windings may result in a negative area. See "Interpreting negative area values" below. + /// - Polygons are assumed to be smaller than half the size of the earth. If you expect to be dealing + /// with polygons larger than this, please use the `unsigned` methods. + /// + /// # Units + /// + /// - return value: meter² + /// + /// # Interpreting negative area values + /// + /// A negative value can mean one of two things: + /// 1. The winding of the polygon is in the clockwise direction (reverse winding). If this is the case, and you know the polygon is smaller than half the area of earth, you can take the absolute value of the reported area to get the correct area. + /// 2. The polygon is larger than half the planet. In this case, the returned area of the polygon is not correct. If you expect to be dealing with very large polygons, please use the `unsigned` methods. + /// + /// [Karney (2013)]: https://arxiv.org/pdf/1109.4448.pdf + #[wasm_bindgen] + pub fn geodesic_area_signed(&self) -> FloatArray { + use geoarrow::algorithm::geo::GeodesicArea; + FloatArray(GeodesicArea::geodesic_area_signed(&self.0)) + } + + /// Determine the area of a geometry on an ellipsoidal model of the earth. Supports very large geometries that cover a significant portion of the earth. + /// + /// This uses the geodesic measurement methods given by [Karney (2013)]. + /// + /// # Assumptions + /// - Polygons are assumed to be wound in a counter-clockwise direction + /// for the exterior ring and a clockwise direction for interior rings. + /// This is the standard winding for geometries that follow the Simple Features standard. + /// Using alternative windings will result in incorrect results. + /// + /// # Units + /// + /// - return value: meter² + /// + /// [Karney (2013)]: https://arxiv.org/pdf/1109.4448.pdf + #[wasm_bindgen] + pub fn geodesic_area_unsigned(&self) -> FloatArray { + use geoarrow::algorithm::geo::GeodesicArea; + FloatArray(GeodesicArea::geodesic_area_unsigned(&self.0)) + } + + /// Determine the perimeter of a geometry on an ellipsoidal model of the earth. + /// + /// This uses the geodesic measurement methods given by [Karney (2013)]. + /// + /// For a polygon this returns the sum of the perimeter of the exterior ring and interior rings. + /// To get the perimeter of just the exterior ring of a polygon, do `polygon.exterior().geodesic_length()`. + /// + /// # Units + /// + /// - return value: meter + /// + /// [Karney (2013)]: https://arxiv.org/pdf/1109.4448.pdf + #[wasm_bindgen] + pub fn geodesic_perimeter(&self) -> FloatArray { + use geoarrow::algorithm::geo::GeodesicArea; + FloatArray(GeodesicArea::geodesic_perimeter(&self.0)) + } + + // TODO: pass tuple of arrays across wasm boundary + + // /// Determine the perimeter and area of a geometry on an ellipsoidal model of the earth, all in one operation. + // /// + // /// This returns the perimeter and area in a `(perimeter, area)` tuple and uses the geodesic measurement methods given by [Karney (2013)]. + // /// + // /// # Area Assumptions + // /// - Polygons are assumed to be wound in a counter-clockwise direction + // /// for the exterior ring and a clockwise direction for interior rings. + // /// This is the standard winding for Geometries that follow the Simple Features standard. + // /// Alternative windings may result in a negative area. See "Interpreting negative area values" below. + // /// - Polygons are assumed to be smaller than half the size of the earth. If you expect to be dealing + // /// with polygons larger than this, please use the 'unsigned' methods. + // /// + // /// # Perimeter + // /// For a polygon this returns the sum of the perimeter of the exterior ring and interior rings. + // /// To get the perimeter of just the exterior ring of a polygon, do `polygon.exterior().geodesic_length()`. + // /// + // /// # Units + // /// + // /// - return value: (meter, meter²) + // /// + // /// # Interpreting negative area values + // /// + // /// A negative area value can mean one of two things: + // /// 1. The winding of the polygon is in the clockwise direction (reverse winding). If this is the case, and you know the polygon is smaller than half the area of earth, you can take the absolute value of the reported area to get the correct area. + // /// 2. The polygon is larger than half the planet. In this case, the returned area of the polygon is not correct. If you expect to be dealing with very large polygons, please use the 'unsigned' methods. + // /// + // /// [Karney (2013)]: https://arxiv.org/pdf/1109.4448.pdf + // #[wasm_bindgen] + // pub fn geodesic_perimeter_area_signed(&self) -> (FloatArray, FloatArray) { + // use geoarrow::algorithm::geo::GeodesicArea; + // let (perimeter, area) = GeodesicArea::geodesic_perimeter_area_signed(&self.0); + // (FloatArray(perimeter), FloatArray(area)) + // } + + // /// Determine the perimeter and area of a geometry on an ellipsoidal model of the earth, all in one operation. Supports very large geometries that cover a significant portion of the earth. + // /// + // /// This returns the perimeter and area in a `(perimeter, area)` tuple and uses the geodesic measurement methods given by [Karney (2013)]. + // /// + // /// # Area Assumptions + // /// - Polygons are assumed to be wound in a counter-clockwise direction + // /// for the exterior ring and a clockwise direction for interior rings. + // /// This is the standard winding for Geometries that follow the Simple Features standard. + // /// Using alternative windings will result in incorrect results. + // /// + // /// # Perimeter + // /// For a polygon this returns the perimeter of the exterior ring and interior rings. + // /// To get the perimeter of just the exterior ring of a polygon, do `polygon.exterior().geodesic_length()`. + // /// + // /// # Units + // /// + // /// - return value: (meter, meter²) + // /// + // /// [Karney (2013)]: https://arxiv.org/pdf/1109.4448.pdf + // fn geodesic_perimeter_area_unsigned(&self) -> (FloatArray, FloatArray); + } + }; +} + +impl_geodesic_area!(PointArray); +impl_geodesic_area!(LineStringArray); +impl_geodesic_area!(PolygonArray); +impl_geodesic_area!(MultiPointArray); +impl_geodesic_area!(MultiLineStringArray); +impl_geodesic_area!(MultiPolygonArray); +impl_geodesic_area!(GeometryArray); diff --git a/js/src/algorithm/geo/mod.rs b/js/src/algorithm/geo/mod.rs index cb7ad4a79..5525e5377 100644 --- a/js/src/algorithm/geo/mod.rs +++ b/js/src/algorithm/geo/mod.rs @@ -5,6 +5,7 @@ pub mod chamberlain_duquette_area; pub mod convex_hull; pub mod dimensions; pub mod euclidean_length; +pub mod geodesic_area; pub mod geodesic_length; pub mod haversine_length; pub mod vincenty_length; diff --git a/js/src/array/geometry.rs b/js/src/array/geometry.rs index 21c207689..6294803a9 100644 --- a/js/src/array/geometry.rs +++ b/js/src/array/geometry.rs @@ -1,6 +1,4 @@ -use crate::array::ffi::FFIArrowArray; use crate::array::polygon::PolygonArray; -use crate::array::primitive::FloatArray; use crate::array::{ LineStringArray, MultiLineStringArray, MultiPointArray, MultiPolygonArray, PointArray, }; @@ -11,8 +9,6 @@ use crate::log; #[cfg(feature = "geodesy")] use crate::reproject::ReprojectDirection; use crate::TransformOrigin; -use arrow2::datatypes::Field; -use geoarrow::GeometryArrayTrait; use wasm_bindgen::prelude::*; #[wasm_bindgen] diff --git a/js/src/array/linestring.rs b/js/src/array/linestring.rs index 6b1dadba4..f8ee291af 100644 --- a/js/src/array/linestring.rs +++ b/js/src/array/linestring.rs @@ -1,6 +1,4 @@ -use crate::array::ffi::FFIArrowArray; use crate::array::primitive::BooleanArray; -use crate::array::primitive::FloatArray; use crate::array::CoordBuffer; use crate::array::GeometryArray; use crate::broadcasting::{BroadcastableAffine, BroadcastableFloat}; @@ -11,8 +9,6 @@ use crate::log; use crate::reproject::ReprojectDirection; use crate::utils::vec_to_offsets; use crate::TransformOrigin; -use arrow2::datatypes::Field; -use geoarrow::GeometryArrayTrait; use wasm_bindgen::prelude::*; #[wasm_bindgen] diff --git a/js/src/array/macro.rs b/js/src/array/macro.rs index b1590e751..0fc68fb49 100644 --- a/js/src/array/macro.rs +++ b/js/src/array/macro.rs @@ -14,18 +14,6 @@ macro_rules! impl_geometry_array { Ok(GeometryArray(affine_transform(&self.into(), transform.0)?)) } - #[wasm_bindgen] - pub fn geodesic_area(&self) -> WasmResult { - use geoarrow::algorithm::geo::geodesic_area_unsigned; - Ok(FloatArray(geodesic_area_unsigned(&self.into())?)) - } - - #[wasm_bindgen] - pub fn geodesic_area_signed(&self) -> WasmResult { - use geoarrow::algorithm::geo::geodesic_area_signed; - Ok(FloatArray(geodesic_area_signed(&self.into())?)) - } - #[cfg(feature = "geodesy")] #[wasm_bindgen] pub fn reproject_rs( @@ -83,13 +71,6 @@ macro_rules! impl_geometry_array { )?)) } - #[wasm_bindgen] - pub fn to_ffi(&self) -> FFIArrowArray { - let arrow_array = self.0.clone().into_boxed_arrow(); - let field = Field::new("", arrow_array.data_type().clone(), true); - FFIArrowArray::new(&field, arrow_array) - } - #[wasm_bindgen] pub fn translate( &self, diff --git a/js/src/array/mod.rs b/js/src/array/mod.rs index d98f2e327..61ddf8c89 100644 --- a/js/src/array/mod.rs +++ b/js/src/array/mod.rs @@ -1,5 +1,4 @@ pub mod coord; -pub mod ffi; pub mod geometry; pub mod linestring; pub mod r#macro; diff --git a/js/src/array/multilinestring.rs b/js/src/array/multilinestring.rs index dc791c865..1e5217b86 100644 --- a/js/src/array/multilinestring.rs +++ b/js/src/array/multilinestring.rs @@ -1,6 +1,4 @@ -use crate::array::ffi::FFIArrowArray; use crate::array::primitive::BooleanArray; -use crate::array::primitive::FloatArray; use crate::array::CoordBuffer; use crate::array::GeometryArray; use crate::broadcasting::{BroadcastableAffine, BroadcastableFloat}; @@ -11,8 +9,6 @@ use crate::log; use crate::reproject::ReprojectDirection; use crate::utils::vec_to_offsets; use crate::TransformOrigin; -use arrow2::datatypes::Field; -use geoarrow::GeometryArrayTrait; use wasm_bindgen::prelude::*; #[wasm_bindgen] diff --git a/js/src/array/multipoint.rs b/js/src/array/multipoint.rs index 7d4a44b55..c8b780459 100644 --- a/js/src/array/multipoint.rs +++ b/js/src/array/multipoint.rs @@ -1,6 +1,4 @@ -use crate::array::ffi::FFIArrowArray; use crate::array::primitive::BooleanArray; -use crate::array::primitive::FloatArray; use crate::array::CoordBuffer; use crate::array::GeometryArray; use crate::broadcasting::{BroadcastableAffine, BroadcastableFloat}; @@ -11,8 +9,6 @@ use crate::log; use crate::reproject::ReprojectDirection; use crate::utils::vec_to_offsets; use crate::TransformOrigin; -use arrow2::datatypes::Field; -use geoarrow::GeometryArrayTrait; use wasm_bindgen::prelude::*; #[wasm_bindgen] diff --git a/js/src/array/multipolygon.rs b/js/src/array/multipolygon.rs index 2debf437a..58840f389 100644 --- a/js/src/array/multipolygon.rs +++ b/js/src/array/multipolygon.rs @@ -1,6 +1,4 @@ -use crate::array::ffi::FFIArrowArray; use crate::array::primitive::BooleanArray; -use crate::array::primitive::FloatArray; use crate::array::CoordBuffer; use crate::array::GeometryArray; use crate::broadcasting::{BroadcastableAffine, BroadcastableFloat}; @@ -11,8 +9,6 @@ use crate::log; use crate::reproject::ReprojectDirection; use crate::utils::vec_to_offsets; use crate::TransformOrigin; -use arrow2::datatypes::Field; -use geoarrow::GeometryArrayTrait; use wasm_bindgen::prelude::*; #[wasm_bindgen] diff --git a/js/src/array/point.rs b/js/src/array/point.rs index cf7dc0f24..571acf73c 100644 --- a/js/src/array/point.rs +++ b/js/src/array/point.rs @@ -1,6 +1,4 @@ -use crate::array::ffi::FFIArrowArray; use crate::array::primitive::BooleanArray; -use crate::array::primitive::FloatArray; use crate::array::CoordBuffer; use crate::array::GeometryArray; use crate::broadcasting::{BroadcastableAffine, BroadcastableFloat}; @@ -10,8 +8,6 @@ use crate::log; #[cfg(feature = "geodesy")] use crate::reproject::ReprojectDirection; use crate::TransformOrigin; -use arrow2::datatypes::Field; -use geoarrow::GeometryArrayTrait; use wasm_bindgen::prelude::*; #[wasm_bindgen] diff --git a/js/src/array/polygon.rs b/js/src/array/polygon.rs index df887b970..f22e68270 100644 --- a/js/src/array/polygon.rs +++ b/js/src/array/polygon.rs @@ -1,6 +1,4 @@ -use crate::array::ffi::FFIArrowArray; use crate::array::primitive::BooleanArray; -use crate::array::primitive::FloatArray; use crate::array::CoordBuffer; use crate::array::GeometryArray; use crate::broadcasting::{BroadcastableAffine, BroadcastableFloat}; @@ -11,8 +9,6 @@ use crate::log; use crate::reproject::ReprojectDirection; use crate::utils::vec_to_offsets; use crate::TransformOrigin; -use arrow2::datatypes::Field; -use geoarrow::GeometryArrayTrait; use wasm_bindgen::prelude::*; #[wasm_bindgen] diff --git a/js/src/array/primitive.rs b/js/src/array/primitive.rs index f0440e606..946ca5288 100644 --- a/js/src/array/primitive.rs +++ b/js/src/array/primitive.rs @@ -1,4 +1,4 @@ -use crate::array::ffi::FFIArrowArray; +use crate::ffi::FFIArrowArray; use arrow2::datatypes::Field; use wasm_bindgen::prelude::*; diff --git a/js/src/array/ffi.rs b/js/src/ffi/array.rs similarity index 100% rename from js/src/array/ffi.rs rename to js/src/ffi/array.rs diff --git a/js/src/ffi/mod.rs b/js/src/ffi/mod.rs new file mode 100644 index 000000000..c81af6955 --- /dev/null +++ b/js/src/ffi/mod.rs @@ -0,0 +1,4 @@ +pub mod array; +pub mod to_ffi; + +pub use array::FFIArrowArray; diff --git a/js/src/ffi/to_ffi.rs b/js/src/ffi/to_ffi.rs new file mode 100644 index 000000000..fae871bc0 --- /dev/null +++ b/js/src/ffi/to_ffi.rs @@ -0,0 +1,27 @@ +use crate::array::*; +use crate::ffi::FFIArrowArray; +use arrow2::datatypes::Field; +use geoarrow::GeometryArrayTrait; +use wasm_bindgen::prelude::*; + +macro_rules! impl_to_ffi { + ($struct_name:ident) => { + #[wasm_bindgen] + impl $struct_name { + #[wasm_bindgen] + pub fn to_ffi(&self) -> FFIArrowArray { + let arrow_array = self.0.clone().into_boxed_arrow(); + let field = Field::new("", arrow_array.data_type().clone(), true); + FFIArrowArray::new(&field, arrow_array) + } + } + }; +} + +impl_to_ffi!(PointArray); +impl_to_ffi!(LineStringArray); +impl_to_ffi!(PolygonArray); +impl_to_ffi!(MultiPointArray); +impl_to_ffi!(MultiLineStringArray); +impl_to_ffi!(MultiPolygonArray); +impl_to_ffi!(GeometryArray); diff --git a/js/src/lib.rs b/js/src/lib.rs index 65f36009f..be6a98fd0 100644 --- a/js/src/lib.rs +++ b/js/src/lib.rs @@ -2,6 +2,7 @@ pub mod algorithm; pub mod array; pub mod broadcasting; pub mod error; +pub mod ffi; #[cfg(feature = "geodesy")] pub mod reproject; pub mod transform_origin; diff --git a/src/algorithm/geo/geodesic_area.rs b/src/algorithm/geo/geodesic_area.rs index ec77270a5..df9579dc8 100644 --- a/src/algorithm/geo/geodesic_area.rs +++ b/src/algorithm/geo/geodesic_area.rs @@ -1,123 +1,295 @@ -use crate::array::GeometryArray; -use crate::error::Result; +use crate::algorithm::geo::utils::zeroes; +use crate::array::*; use crate::GeometryArrayTrait; use arrow2::array::{MutablePrimitiveArray, PrimitiveArray}; -use geo::prelude::GeodesicArea; +use geo::prelude::GeodesicArea as _GeodesicArea; -/// Calculate the unsigned geodesic area of a geometry on an ellipsoid using the algorithm -/// presented in Algorithms for geodesics by Charles Karney (2013) -pub fn geodesic_area_unsigned(array: &GeometryArray) -> Result> { - let mut output_array = MutablePrimitiveArray::::with_capacity(array.len()); +/// Determine the perimeter and area of a geometry on an ellipsoidal model of the earth. +/// +/// This uses the geodesic measurement methods given by [Karney (2013)]. +/// +/// [Karney (2013)]: https://arxiv.org/pdf/1109.4448.pdf +pub trait GeodesicArea { + /// Determine the area of a geometry on an ellipsoidal model of the earth. + /// + /// This uses the geodesic measurement methods given by [Karney (2013)]. + /// + /// # Assumptions + /// - Polygons are assumed to be wound in a counter-clockwise direction + /// for the exterior ring and a clockwise direction for interior rings. + /// This is the standard winding for geometries that follow the Simple Feature standard. + /// Alternative windings may result in a negative area. See "Interpreting negative area values" below. + /// - Polygons are assumed to be smaller than half the size of the earth. If you expect to be dealing + /// with polygons larger than this, please use the `unsigned` methods. + /// + /// # Units + /// + /// - return value: meter² + /// + /// # Interpreting negative area values + /// + /// A negative value can mean one of two things: + /// 1. The winding of the polygon is in the clockwise direction (reverse winding). If this is the case, and you know the polygon is smaller than half the area of earth, you can take the absolute value of the reported area to get the correct area. + /// 2. The polygon is larger than half the planet. In this case, the returned area of the polygon is not correct. If you expect to be dealing with very large polygons, please use the `unsigned` methods. + /// + /// # Examples + /// ```rust + /// use geo::prelude::*; + /// use geo::polygon; + /// use geo::Polygon; + /// + /// // The O2 in London + /// let mut polygon: Polygon = polygon![ + /// (x: 0.00388383, y: 51.501574), + /// (x: 0.00538587, y: 51.502278), + /// (x: 0.00553607, y: 51.503299), + /// (x: 0.00467777, y: 51.504181), + /// (x: 0.00327229, y: 51.504435), + /// (x: 0.00187754, y: 51.504168), + /// (x: 0.00087976, y: 51.503380), + /// (x: 0.00107288, y: 51.502324), + /// (x: 0.00185608, y: 51.501770), + /// (x: 0.00388383, y: 51.501574), + /// ]; + /// + /// let area = polygon.geodesic_area_unsigned(); + /// + /// assert_eq!( + /// 78_596., // meters + /// area.round() + /// ); + /// ``` + /// [Karney (2013)]: https://arxiv.org/pdf/1109.4448.pdf + fn geodesic_area_signed(&self) -> PrimitiveArray; - match array { - GeometryArray::WKB(arr) => { - arr.iter_geo() - .for_each(|maybe_g| output_array.push(maybe_g.map(|g| g.geodesic_area_unsigned()))); - } - GeometryArray::Point(arr) => { - arr.iter_geo() - .for_each(|maybe_g| output_array.push(maybe_g.map(|g| g.geodesic_area_unsigned()))); - } - GeometryArray::LineString(arr) => { - arr.iter_geo() - .for_each(|maybe_g| output_array.push(maybe_g.map(|g| g.geodesic_area_unsigned()))); - } - GeometryArray::Polygon(arr) => { - arr.iter_geo() - .for_each(|maybe_g| output_array.push(maybe_g.map(|g| g.geodesic_area_unsigned()))); - } - GeometryArray::MultiPoint(arr) => { - arr.iter_geo() - .for_each(|maybe_g| output_array.push(maybe_g.map(|g| g.geodesic_area_unsigned()))); - } - GeometryArray::MultiLineString(arr) => { - arr.iter_geo() - .for_each(|maybe_g| output_array.push(maybe_g.map(|g| g.geodesic_area_unsigned()))); - } - GeometryArray::MultiPolygon(arr) => { - arr.iter_geo() - .for_each(|maybe_g| output_array.push(maybe_g.map(|g| g.geodesic_area_unsigned()))); - } - } + /// Determine the area of a geometry on an ellipsoidal model of the earth. Supports very large geometries that cover a significant portion of the earth. + /// + /// This uses the geodesic measurement methods given by [Karney (2013)]. + /// + /// # Assumptions + /// - Polygons are assumed to be wound in a counter-clockwise direction + /// for the exterior ring and a clockwise direction for interior rings. + /// This is the standard winding for geometries that follow the Simple Features standard. + /// Using alternative windings will result in incorrect results. + /// + /// # Units + /// + /// - return value: meter² + /// + /// # Examples + /// ```rust + /// use geo::prelude::*; + /// use geo::polygon; + /// use geo::Polygon; + /// + /// // Describe a polygon that covers all of the earth EXCEPT this small square. + /// // The outside of the polygon is in this square, the inside of the polygon is the rest of the earth. + /// let mut polygon: Polygon = polygon![ + /// (x: 0.0, y: 0.0), + /// (x: 0.0, y: 1.0), + /// (x: 1.0, y: 1.0), + /// (x: 1.0, y: 0.0), + /// ]; + /// + /// let area = polygon.geodesic_area_unsigned(); + /// + /// // Over 5 trillion square meters! + /// assert_eq!(area, 510053312945726.94); + /// ``` + /// [Karney (2013)]: https://arxiv.org/pdf/1109.4448.pdf + fn geodesic_area_unsigned(&self) -> PrimitiveArray; - Ok(output_array.into()) + /// Determine the perimeter of a geometry on an ellipsoidal model of the earth. + /// + /// This uses the geodesic measurement methods given by [Karney (2013)]. + /// + /// For a polygon this returns the sum of the perimeter of the exterior ring and interior rings. + /// To get the perimeter of just the exterior ring of a polygon, do `polygon.exterior().geodesic_length()`. + /// + /// # Units + /// + /// - return value: meter + /// + /// [Karney (2013)]: https://arxiv.org/pdf/1109.4448.pdf + fn geodesic_perimeter(&self) -> PrimitiveArray; + + /// Determine the perimeter and area of a geometry on an ellipsoidal model of the earth, all in one operation. + /// + /// This returns the perimeter and area in a `(perimeter, area)` tuple and uses the geodesic measurement methods given by [Karney (2013)]. + /// + /// # Area Assumptions + /// - Polygons are assumed to be wound in a counter-clockwise direction + /// for the exterior ring and a clockwise direction for interior rings. + /// This is the standard winding for Geometries that follow the Simple Features standard. + /// Alternative windings may result in a negative area. See "Interpreting negative area values" below. + /// - Polygons are assumed to be smaller than half the size of the earth. If you expect to be dealing + /// with polygons larger than this, please use the 'unsigned' methods. + /// + /// # Perimeter + /// For a polygon this returns the sum of the perimeter of the exterior ring and interior rings. + /// To get the perimeter of just the exterior ring of a polygon, do `polygon.exterior().geodesic_length()`. + /// + /// # Units + /// + /// - return value: (meter, meter²) + /// + /// # Interpreting negative area values + /// + /// A negative area value can mean one of two things: + /// 1. The winding of the polygon is in the clockwise direction (reverse winding). If this is the case, and you know the polygon is smaller than half the area of earth, you can take the absolute value of the reported area to get the correct area. + /// 2. The polygon is larger than half the planet. In this case, the returned area of the polygon is not correct. If you expect to be dealing with very large polygons, please use the 'unsigned' methods. + /// + /// [Karney (2013)]: https://arxiv.org/pdf/1109.4448.pdf + fn geodesic_perimeter_area_signed(&self) -> (PrimitiveArray, PrimitiveArray); + + /// Determine the perimeter and area of a geometry on an ellipsoidal model of the earth, all in one operation. Supports very large geometries that cover a significant portion of the earth. + /// + /// This returns the perimeter and area in a `(perimeter, area)` tuple and uses the geodesic measurement methods given by [Karney (2013)]. + /// + /// # Area Assumptions + /// - Polygons are assumed to be wound in a counter-clockwise direction + /// for the exterior ring and a clockwise direction for interior rings. + /// This is the standard winding for Geometries that follow the Simple Features standard. + /// Using alternative windings will result in incorrect results. + /// + /// # Perimeter + /// For a polygon this returns the perimeter of the exterior ring and interior rings. + /// To get the perimeter of just the exterior ring of a polygon, do `polygon.exterior().geodesic_length()`. + /// + /// # Units + /// + /// - return value: (meter, meter²) + /// + /// [Karney (2013)]: https://arxiv.org/pdf/1109.4448.pdf + fn geodesic_perimeter_area_unsigned(&self) -> (PrimitiveArray, PrimitiveArray); } -/// Calculate the signed geodesic area of a geometry on an ellipsoid using the algorithm -/// presented in Algorithms for geodesics by Charles Karney (2013) -pub fn geodesic_area_signed(array: &GeometryArray) -> Result> { - let mut output_array = MutablePrimitiveArray::::with_capacity(array.len()); +/// Generate a `GeodesicArea` implementation where the result is zero. +macro_rules! zero_impl { + ($type:ident) => { + impl GeodesicArea for $type { + fn geodesic_perimeter(&self) -> PrimitiveArray { + zeroes(self.len(), self.validity()) + } - match array { - GeometryArray::WKB(arr) => { - arr.iter_geo() - .for_each(|maybe_g| output_array.push(maybe_g.map(|g| g.geodesic_area_signed()))); - } - GeometryArray::Point(arr) => { - arr.iter_geo() - .for_each(|maybe_g| output_array.push(maybe_g.map(|g| g.geodesic_area_signed()))); - } - GeometryArray::LineString(arr) => { - arr.iter_geo() - .for_each(|maybe_g| output_array.push(maybe_g.map(|g| g.geodesic_area_signed()))); - } - GeometryArray::Polygon(arr) => { - arr.iter_geo() - .for_each(|maybe_g| output_array.push(maybe_g.map(|g| g.geodesic_area_signed()))); - } - GeometryArray::MultiPoint(arr) => { - arr.iter_geo() - .for_each(|maybe_g| output_array.push(maybe_g.map(|g| g.geodesic_area_signed()))); - } - GeometryArray::MultiLineString(arr) => { - arr.iter_geo() - .for_each(|maybe_g| output_array.push(maybe_g.map(|g| g.geodesic_area_signed()))); - } - GeometryArray::MultiPolygon(arr) => { - arr.iter_geo() - .for_each(|maybe_g| output_array.push(maybe_g.map(|g| g.geodesic_area_signed()))); - } - } + fn geodesic_area_signed(&self) -> PrimitiveArray { + zeroes(self.len(), self.validity()) + } - Ok(output_array.into()) -} + fn geodesic_area_unsigned(&self) -> PrimitiveArray { + zeroes(self.len(), self.validity()) + } -/// Determine the perimeter of a geometry on an ellipsoidal model of the earth. -/// -/// This uses the geodesic measurement methods given by Karney (2013). -pub fn geodesic_perimeter(array: GeometryArray) -> Result> { - let mut output_array = MutablePrimitiveArray::::with_capacity(array.len()); - - match array { - GeometryArray::WKB(arr) => { - arr.iter_geo() - .for_each(|maybe_g| output_array.push(maybe_g.map(|g| g.geodesic_perimeter()))); - } - GeometryArray::Point(arr) => { - arr.iter_geo() - .for_each(|maybe_g| output_array.push(maybe_g.map(|g| g.geodesic_perimeter()))); - } - GeometryArray::LineString(arr) => { - arr.iter_geo() - .for_each(|maybe_g| output_array.push(maybe_g.map(|g| g.geodesic_perimeter()))); - } - GeometryArray::Polygon(arr) => { - arr.iter_geo() - .for_each(|maybe_g| output_array.push(maybe_g.map(|g| g.geodesic_perimeter()))); - } - GeometryArray::MultiPoint(arr) => { - arr.iter_geo() - .for_each(|maybe_g| output_array.push(maybe_g.map(|g| g.geodesic_perimeter()))); - } - GeometryArray::MultiLineString(arr) => { - arr.iter_geo() - .for_each(|maybe_g| output_array.push(maybe_g.map(|g| g.geodesic_perimeter()))); + fn geodesic_perimeter_area_signed(&self) -> (PrimitiveArray, PrimitiveArray) { + ( + zeroes(self.len(), self.validity()), + zeroes(self.len(), self.validity()), + ) + } + + fn geodesic_perimeter_area_unsigned( + &self, + ) -> (PrimitiveArray, PrimitiveArray) { + ( + zeroes(self.len(), self.validity()), + zeroes(self.len(), self.validity()), + ) + } } - GeometryArray::MultiPolygon(arr) => { - arr.iter_geo() - .for_each(|maybe_g| output_array.push(maybe_g.map(|g| g.geodesic_perimeter()))); + }; +} + +zero_impl!(PointArray); +zero_impl!(LineStringArray); +zero_impl!(MultiPointArray); +zero_impl!(MultiLineStringArray); + +/// Implementation that iterates over geo objects +macro_rules! iter_geo_impl { + ($type:ident) => { + impl GeodesicArea for $type { + fn geodesic_perimeter(&self) -> PrimitiveArray { + let mut output_array = MutablePrimitiveArray::::with_capacity(self.len()); + + self.iter_geo() + .for_each(|maybe_g| output_array.push(maybe_g.map(|g| g.geodesic_perimeter()))); + + output_array.into() + } + + fn geodesic_area_signed(&self) -> PrimitiveArray { + let mut output_array = MutablePrimitiveArray::::with_capacity(self.len()); + + self.iter_geo().for_each(|maybe_g| { + output_array.push(maybe_g.map(|g| g.geodesic_area_signed())) + }); + + output_array.into() + } + + fn geodesic_area_unsigned(&self) -> PrimitiveArray { + let mut output_array = MutablePrimitiveArray::::with_capacity(self.len()); + + self.iter_geo().for_each(|maybe_g| { + output_array.push(maybe_g.map(|g| g.geodesic_area_unsigned())) + }); + + output_array.into() + } + + fn geodesic_perimeter_area_signed(&self) -> (PrimitiveArray, PrimitiveArray) { + let mut output_perimeter_array = + MutablePrimitiveArray::::with_capacity(self.len()); + let mut output_area_array = MutablePrimitiveArray::::with_capacity(self.len()); + + self.iter_geo().for_each(|maybe_g| { + if let Some(g) = maybe_g { + let (perimeter, area) = g.geodesic_perimeter_area_signed(); + output_perimeter_array.push(Some(perimeter)); + output_area_array.push(Some(area)); + } else { + output_perimeter_array.push(None); + output_area_array.push(None); + } + }); + + (output_perimeter_array.into(), output_area_array.into()) + } + + fn geodesic_perimeter_area_unsigned( + &self, + ) -> (PrimitiveArray, PrimitiveArray) { + let mut output_perimeter_array = + MutablePrimitiveArray::::with_capacity(self.len()); + let mut output_area_array = MutablePrimitiveArray::::with_capacity(self.len()); + + self.iter_geo().for_each(|maybe_g| { + if let Some(g) = maybe_g { + let (perimeter, area) = g.geodesic_perimeter_area_unsigned(); + output_perimeter_array.push(Some(perimeter)); + output_area_array.push(Some(area)); + } else { + output_perimeter_array.push(None); + output_area_array.push(None); + } + }); + + (output_perimeter_array.into(), output_area_array.into()) + } } - } + }; +} + +iter_geo_impl!(PolygonArray); +iter_geo_impl!(MultiPolygonArray); +iter_geo_impl!(WKBArray); - Ok(output_array.into()) +impl GeodesicArea for GeometryArray { + crate::geometry_array_delegate_impl! { + fn geodesic_area_signed(&self) -> PrimitiveArray; + fn geodesic_area_unsigned(&self) -> PrimitiveArray; + fn geodesic_perimeter(&self) -> PrimitiveArray; + fn geodesic_perimeter_area_signed(&self) -> (PrimitiveArray, PrimitiveArray); + fn geodesic_perimeter_area_unsigned(&self) -> (PrimitiveArray, PrimitiveArray); + } } diff --git a/src/algorithm/geo/mod.rs b/src/algorithm/geo/mod.rs index fef42cdfd..ac0b4dabc 100644 --- a/src/algorithm/geo/mod.rs +++ b/src/algorithm/geo/mod.rs @@ -2,12 +2,10 @@ mod affine; mod distance; -mod geodesic_area; mod simplify; pub(crate) mod utils; pub use affine::{affine_transform, rotate, scale, skew, translate, TransformOrigin}; -pub use geodesic_area::{geodesic_area_signed, geodesic_area_unsigned, geodesic_perimeter}; pub use simplify::simplify; /// Composable affine operations such as rotate, scale, skew, and translate @@ -47,6 +45,10 @@ pub use dimensions::HasDimensions; pub mod euclidean_length; pub use euclidean_length::EuclideanLength; +/// Calculate the Geodesic area and perimeter of polygons. +pub mod geodesic_area; +pub use geodesic_area::GeodesicArea; + /// Calculate the Geodesic length of a line. pub mod geodesic_length; pub use geodesic_length::GeodesicLength;