From 652df618b57aec34e8d598ae69d0830376fda5b3 Mon Sep 17 00:00:00 2001 From: Grzegorz Zakrzewski Date: Thu, 18 Sep 2025 20:34:26 +0000 Subject: [PATCH] feat: oblique lon lat transformation --- src/projections.rs | 2 + src/projections/oblique_lon_lat.rs | 127 +++++++++++++++++++++++++ tests/basic_correctness.rs | 5 + tests/special_cases/mod.rs | 1 + tests/special_cases/oblique_lon_lat.rs | 67 +++++++++++++ 5 files changed, 202 insertions(+) create mode 100644 src/projections/oblique_lon_lat.rs create mode 100644 tests/special_cases/oblique_lon_lat.rs diff --git a/src/projections.rs b/src/projections.rs index 98144ef..dfae0dd 100644 --- a/src/projections.rs +++ b/src/projections.rs @@ -5,9 +5,11 @@ mod lambert_conformal_conic; mod modified_azimuthal_equidistant; mod equidistant_cylindrical; mod lon_lat; +mod oblique_lon_lat; pub use azimuthal_equidistant::AzimuthalEquidistant; pub use lambert_conformal_conic::LambertConformalConic; pub use modified_azimuthal_equidistant::ModifiedAzimuthalEquidistant; pub use equidistant_cylindrical::EquidistantCylindrical; pub use lon_lat::LongitudeLatitude; +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..8da6a1c --- /dev/null +++ b/src/projections/oblique_lon_lat.rs @@ -0,0 +1,127 @@ +use crate::Projection; +use crate::errors::ProjectionError; + +/// An oblique longitude-latitude projection 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. +#[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 { + /// 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. + /// + /// # Arguments + /// - `pole_lon`, `pole_lat` - longitude and latitude of the original North Pole (90°N, 0°E) in the new rotated system + /// - `central_lon` - longitude of the central meridian (optional, defaults to 0.0) + /// + /// # Errors + /// + /// Returns [`ProjectionError::IncorrectParams`] 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 new(pole_lon: f64, pole_lat: f64, central_lon: Option) -> Result { + let central_lon = central_lon.unwrap_or(0.0); + if !pole_lon.is_finite() || !pole_lat.is_finite() || !central_lon.is_finite() { + return Err(ProjectionError::IncorrectParams( + "one of arguments is not finite", + )); + } + + if !(-180.0..180.0).contains(&pole_lon) || !(-180.0..180.0).contains(¢ral_lon) { + return Err(ProjectionError::IncorrectParams( + "longitude must be between -180..180", + )); + } + + if !(-90.0..90.0).contains(&pole_lat) { + return Err(ProjectionError::IncorrectParams( + "latitude must be between -90..90", + )); + } + + 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 { + 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); + } + + 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 c48718b..bacc58d 100644 --- a/tests/basic_correctness.rs +++ b/tests/basic_correctness.rs @@ -127,6 +127,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..0cef0b6 --- /dev/null +++ b/tests/special_cases/oblique_lon_lat.rs @@ -0,0 +1,67 @@ +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::new(pole_lon as f64, pole_lat as f64, Some(central_lon as f64)) + .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); + } +}