diff --git a/src/projections.rs b/src/projections.rs index 0815f3c..2b55633 100644 --- a/src/projections.rs +++ b/src/projections.rs @@ -4,6 +4,7 @@ pub mod azimuthal_equidistant; pub mod equidistant_cylindrical; pub mod lambert_conformal_conic; mod lon_lat; +mod oblique_lon_lat; pub mod modified_azimuthal_equidistant; pub use azimuthal_equidistant::AzimuthalEquidistant; @@ -11,3 +12,4 @@ pub use equidistant_cylindrical::EquidistantCylindrical; pub use lambert_conformal_conic::LambertConformalConic; pub use lon_lat::LongitudeLatitude; pub use modified_azimuthal_equidistant::ModifiedAzimuthalEquidistant; +pub use oblique_lon_lat::ObliqueLonLat; diff --git a/src/projections/oblique_lon_lat.rs b/src/projections/oblique_lon_lat.rs new file mode 100644 index 0000000..e19518e --- /dev/null +++ b/src/projections/oblique_lon_lat.rs @@ -0,0 +1,171 @@ +//! An oblique longitude-latitude projection (also known as Rotated Pole) that +//! transforms the Earth's graticule so that the "north pole" of the coordinate system +//! assumes a position different from the true North Pole. This effectively rotates the +//! meridian and parallel grid, allowing any point on Earth to serve as the reference +//! pole of the projection. +//! +//! The transformation is applied by rotating the spherical coordinate system before +//! applying standard longitude-latitude coordinates. This is particularly useful for +//! regional mapping where the area of interest can be positioned optimally relative +//! to the coordinate grid, reducing distortion in the region of interest. +//! +//! Unlike most other projections in this crate, the output of the +//! [`project`](crate::Projection::project) function is in degrees, not meters. + +use crate::Projection; +use crate::errors::ProjectionError; +use crate::errors::{ensure_finite, ensure_within_range, unpack_required_parameter}; + +#[cfg(feature = "tracing")] +use tracing::instrument; + +/// Main projection struct that is constructed from [`ObliqueLonLatBuilder`] and used for computations. +#[derive(Copy, Clone, PartialEq, PartialOrd, Debug)] +pub struct ObliqueLonLat { + lambda_p: f64, + sin_phi_p: f64, + cos_phi_p: f64, + lon_0: f64, +} + +impl ObliqueLonLat { + /// Initializes builder with default values. + /// Projection parameters can be set with builder methods, + /// refer to the documentation of those methods to check which parameters are required + /// and default values for optional arguments. + #[must_use] + pub fn builder() -> ObliqueLonLatBuilder { + ObliqueLonLatBuilder::default() + } +} + +/// Builder struct which allows to construct [`ObliqueLonLat`] projection. +/// Refer to the documentation of this struct's methods to check which parameters are required +/// and default values for optional arguments. +#[derive(Debug)] +pub struct ObliqueLonLatBuilder { + pole_lon: Option, + pole_lat: Option, + central_lon: f64, +} + +impl Default for ObliqueLonLatBuilder { + fn default() -> Self { + Self { + pole_lon: None, + pole_lat: None, + central_lon: 0.0, + } + } +} + +impl ObliqueLonLatBuilder { + /// *(required)* Sets the longitude and latitude of the rotated pole. + pub fn pole_lonlat(&mut self, lon: f64, lat: f64) -> &mut Self { + self.pole_lon = Some(lon); + self.pole_lat = Some(lat); + self + } + + /// *(optional)* Sets the central meridian longitude, defaults to `0.0`. + pub fn central_lon(&mut self, lon: f64) -> &mut Self { + self.central_lon = lon; + self + } + + /// ObliqueLonLat projection constructor. + /// + /// To reduce computational overhead of projection functions this + /// constructor is non-trivial and tries to do as much projection computations as possible. + /// Thus creating a new structure can involve a significant computational overhead. + /// When projecting multiple coordinates only one instance of the structure should be created + /// and copied/borrowed as needed. + /// + /// Ellipsoid is not definable as this projection does not depend on it. + /// + /// # Errors + /// + /// Returns [`ProjectionError`] with additional information when: + /// + /// - one or more longitudes are not within -180..180 range. + /// - one or more latitudes are not within -90..90 range. + /// - one or more arguments are not finite. + pub fn initialize_projection(&self) -> Result { + let pole_lon = unpack_required_parameter!(self, pole_lon); + let pole_lat = unpack_required_parameter!(self, pole_lat); + let central_lon = self.central_lon; + + ensure_finite!(pole_lon, pole_lat, central_lon); + ensure_within_range!(pole_lon, -180.0..180.0); + ensure_within_range!(pole_lat, -90.0..90.0); + ensure_within_range!(central_lon, -180.0..180.0); + + let phi_p = pole_lat.to_radians(); + + Ok(ObliqueLonLat { + lambda_p: pole_lon.to_radians(), + sin_phi_p: phi_p.sin(), + cos_phi_p: phi_p.cos(), + lon_0: central_lon, + }) + } +} + +fn adjust_lon(lon: f64) -> f64 { + let pi_degrees = 180.0_f64; + if lon > pi_degrees { + lon - 2.0 * pi_degrees + } else if lon < -pi_degrees { + lon + 2.0 * pi_degrees + } else { + lon + } +} + +impl Projection for ObliqueLonLat { + #[inline] + #[cfg_attr(feature = "tracing", instrument(level = "trace"))] + fn project_unchecked(&self, lon: f64, lat: f64) -> (f64, f64) { + let lambda = (lon - self.lon_0).to_radians(); + let phi = lat.to_radians(); + + let cos_lambda = lambda.cos(); + let sin_lambda = lambda.sin(); + let cos_phi = phi.cos(); + let sin_phi = phi.sin(); + + // Formula (5-8b) of [Snyder (1987)](https://pubs.er.usgs.gov/publication/pp1395) + let lambda_prime = (cos_phi * sin_lambda) + .atan2(self.sin_phi_p * cos_phi * cos_lambda + self.cos_phi_p * sin_phi) + self.lambda_p; + + // Formula (5-7) + let phi_prime = (self.sin_phi_p * sin_phi - self.cos_phi_p * cos_phi * cos_lambda).asin(); + + let lon_prime = adjust_lon(lambda_prime.to_degrees()); + let lat_prime = phi_prime.to_degrees(); + return (lon_prime, lat_prime); + } + + #[inline] + #[cfg_attr(feature = "tracing", instrument(level = "trace"))] + fn inverse_project_unchecked(&self, x: f64, y: f64) -> (f64, f64) { + let lambda_prime = x.to_radians() - self.lambda_p; + let phi_prime = y.to_radians(); + + let cos_lambda_prime = lambda_prime.cos(); + let sin_lambda_prime = lambda_prime.sin(); + let cos_phi_prime = phi_prime.cos(); + let sin_phi_prime = phi_prime.sin(); + + // Formula (5-10b) + let lambda = (cos_phi_prime * sin_lambda_prime) + .atan2(self.sin_phi_p * cos_phi_prime * cos_lambda_prime - self.cos_phi_p * sin_phi_prime); + + // Formula (5-9) + let phi = (self.sin_phi_p * sin_phi_prime + self.cos_phi_p * cos_phi_prime * cos_lambda_prime).asin(); + + let lon = adjust_lon(lambda.to_degrees() + self.lon_0); + let lat = phi.to_degrees(); + return (lon, lat); + } +} diff --git a/tests/basic_correctness.rs b/tests/basic_correctness.rs index d278b5a..4e93ae1 100644 --- a/tests/basic_correctness.rs +++ b/tests/basic_correctness.rs @@ -114,6 +114,11 @@ fn equidistant_cylindrical() { special_cases::equidistant_cylindrical::basic_correctness(); } +#[test] +fn oblique_lon_lat() { + special_cases::oblique_lon_lat::basic_correctness(); +} + pub fn test_points_with_proj(int_proj: &P, proj_str: &str, extent: TestExtent) { let ref_proj = Proj::new(proj_str).unwrap(); diff --git a/tests/special_cases/mod.rs b/tests/special_cases/mod.rs index cec4aa1..f07d8b2 100644 --- a/tests/special_cases/mod.rs +++ b/tests/special_cases/mod.rs @@ -1,3 +1,4 @@ pub(crate) mod equidistant_cylindrical; pub(crate) mod modified_azimuthal_equidistant; pub(crate) mod lambert_conformal_conic; +pub(crate) mod oblique_lon_lat; diff --git a/tests/special_cases/oblique_lon_lat.rs b/tests/special_cases/oblique_lon_lat.rs new file mode 100644 index 0000000..524a875 --- /dev/null +++ b/tests/special_cases/oblique_lon_lat.rs @@ -0,0 +1,69 @@ +use float_cmp::assert_approx_eq; +use mappers::projections::ObliqueLonLat; +use mappers::Projection; +use proj::Proj; +use crate::GLOBAL_GEO_POINTS; +use crate::LOCAL_GEO_POINTS; +use crate::TestExtent; + +pub(crate) fn basic_correctness() { + // This projection does not depend on ellipsoid + // Also, this projection produces degrees, and not meters, so it must be tested separately + + for pole_lon in (-180..180).step_by(60) { + for pole_lat in (-90..90).step_by(60) { + for central_lon in (-180..180).step_by(60) { + let int_proj = ObliqueLonLat::builder() + .pole_lonlat(pole_lon as f64, pole_lat as f64) + .central_lon(central_lon as f64) + .initialize_projection() + .unwrap(); + + let proj_str = format!( + "+ellps=sphere +proj=ob_tran +o_proj=latlon +o_lat_p={} +o_lon_p={} +lon_0={}", + pole_lat, pole_lon, central_lon + ); + + test_points_with_proj(&int_proj, &proj_str, TestExtent::Global); + test_points_with_proj(&int_proj, &proj_str, TestExtent::Local); + } + } + } +} + +fn test_points_with_proj(int_proj: &ObliqueLonLat, proj_str: &str, extent: TestExtent) { + let ref_proj = Proj::new(proj_str).unwrap(); + + let geo_points = match extent { + TestExtent::Global => GLOBAL_GEO_POINTS, + TestExtent::Local => LOCAL_GEO_POINTS, + }; + + for point in geo_points { + // projection + let (ref_x, ref_y) = ref_proj + .project((point.0.to_radians(), point.1.to_radians()), false) + .unwrap(); + + let ref_lon = ref_x.to_degrees(); + let ref_lat = ref_y.to_degrees(); + + let (tst_lon, tst_lat) = int_proj.project(point.0, point.1).unwrap(); + + assert_approx_eq!(f64, ref_lon, tst_lon, epsilon = 0.000_000_1); + assert_approx_eq!(f64, ref_lat, tst_lat, epsilon = 0.000_000_1); + + // inverse projection + let (ref_x, ref_y) = ref_proj + .project((point.0.to_radians(), point.1.to_radians()), true) + .unwrap(); + + let ref_lon = ref_x.to_degrees(); + let ref_lat = ref_y.to_degrees(); + + let (tst_lon, tst_lat) = int_proj.inverse_project(point.0, point.1).unwrap(); + + assert_approx_eq!(f64, ref_lon, tst_lon, epsilon = 0.000_000_1); + assert_approx_eq!(f64, ref_lat, tst_lat, epsilon = 0.000_000_1); + } +}