Skip to main content

proj_core/
crs.rs

1use crate::datum::Datum;
2use crate::error::{Error, Result};
3
4/// A coordinate system's projected linear unit.
5///
6/// The stored value is the conversion factor from one native unit to meters.
7#[derive(Debug, Clone, Copy, PartialEq)]
8pub struct LinearUnit {
9    meters_per_unit: f64,
10}
11
12impl LinearUnit {
13    /// Metre-based projected coordinates.
14    pub const fn metre() -> Self {
15        Self {
16            meters_per_unit: 1.0,
17        }
18    }
19
20    /// Alias for [`LinearUnit::metre`].
21    pub const fn meter() -> Self {
22        Self::metre()
23    }
24
25    /// Kilometer-based projected coordinates.
26    pub const fn kilometre() -> Self {
27        Self {
28            meters_per_unit: 1000.0,
29        }
30    }
31
32    /// Alias for [`LinearUnit::kilometre`].
33    pub const fn kilometer() -> Self {
34        Self::kilometre()
35    }
36
37    /// International foot-based projected coordinates.
38    pub const fn foot() -> Self {
39        Self {
40            meters_per_unit: 0.3048,
41        }
42    }
43
44    /// US survey foot-based projected coordinates.
45    pub const fn us_survey_foot() -> Self {
46        Self {
47            meters_per_unit: 0.3048006096012192,
48        }
49    }
50
51    /// Construct a custom projected linear unit from its meter conversion factor.
52    pub fn from_meters_per_unit(meters_per_unit: f64) -> Result<Self> {
53        if !meters_per_unit.is_finite() || meters_per_unit <= 0.0 {
54            return Err(Error::InvalidDefinition(
55                "linear unit conversion factor must be a finite positive number".into(),
56            ));
57        }
58
59        Ok(Self { meters_per_unit })
60    }
61
62    /// Return the number of meters represented by one native projected unit.
63    pub const fn meters_per_unit(self) -> f64 {
64        self.meters_per_unit
65    }
66
67    /// Convert a native projected coordinate value into meters.
68    pub const fn to_meters(self, value: f64) -> f64 {
69        value * self.meters_per_unit
70    }
71
72    /// Convert a meter value into the native projected unit.
73    pub const fn from_meters(self, value: f64) -> f64 {
74        value / self.meters_per_unit
75    }
76}
77
78/// A Coordinate Reference System definition.
79#[derive(Debug, Clone)]
80pub enum CrsDef {
81    /// Geographic CRS (lon/lat in degrees).
82    Geographic(GeographicCrsDef),
83    /// Projected CRS (easting/northing in the CRS's native linear unit).
84    Projected(ProjectedCrsDef),
85}
86
87impl CrsDef {
88    /// Get the datum for this CRS.
89    pub fn datum(&self) -> &Datum {
90        match self {
91            CrsDef::Geographic(g) => g.datum(),
92            CrsDef::Projected(p) => p.datum(),
93        }
94    }
95
96    /// Get the EPSG code for this CRS.
97    pub fn epsg(&self) -> u32 {
98        match self {
99            CrsDef::Geographic(g) => g.epsg(),
100            CrsDef::Projected(p) => p.epsg(),
101        }
102    }
103
104    /// Get the CRS name.
105    pub fn name(&self) -> &str {
106        match self {
107            CrsDef::Geographic(g) => g.name(),
108            CrsDef::Projected(p) => p.name(),
109        }
110    }
111
112    /// Returns true if this is a geographic CRS.
113    pub fn is_geographic(&self) -> bool {
114        matches!(self, CrsDef::Geographic(_))
115    }
116
117    /// Returns true if this is a projected CRS.
118    pub fn is_projected(&self) -> bool {
119        matches!(self, CrsDef::Projected(_))
120    }
121
122    /// Returns the geographic CRS EPSG code used for operation selection, when known.
123    pub fn base_geographic_crs_epsg(&self) -> Option<u32> {
124        match self {
125            CrsDef::Geographic(g) if g.epsg() != 0 => Some(g.epsg()),
126            CrsDef::Projected(p) if p.base_geographic_crs_epsg() != 0 => {
127                Some(p.base_geographic_crs_epsg())
128            }
129            _ => None,
130        }
131    }
132
133    /// Returns true when two CRS definitions map to the same internal semantics.
134    pub fn semantically_equivalent(&self, other: &Self) -> bool {
135        match (self, other) {
136            (CrsDef::Geographic(a), CrsDef::Geographic(b)) => a.datum().same_datum(b.datum()),
137            (CrsDef::Projected(a), CrsDef::Projected(b)) => {
138                a.datum().same_datum(b.datum())
139                    && approx_eq(a.linear_unit_to_meter(), b.linear_unit_to_meter())
140                    && projection_methods_equivalent(&a.method(), &b.method())
141            }
142            _ => false,
143        }
144    }
145}
146
147/// Definition of a geographic CRS (longitude, latitude in degrees).
148#[derive(Debug, Clone)]
149pub struct GeographicCrsDef {
150    epsg: u32,
151    datum: Datum,
152    name: &'static str,
153}
154
155impl GeographicCrsDef {
156    pub const fn new(epsg: u32, datum: Datum, name: &'static str) -> Self {
157        Self { epsg, datum, name }
158    }
159
160    pub const fn epsg(&self) -> u32 {
161        self.epsg
162    }
163
164    pub const fn datum(&self) -> &Datum {
165        &self.datum
166    }
167
168    pub const fn name(&self) -> &'static str {
169        self.name
170    }
171}
172
173/// Definition of a projected CRS (easting, northing in the CRS's native linear unit).
174#[derive(Debug, Clone)]
175pub struct ProjectedCrsDef {
176    epsg: u32,
177    base_geographic_crs_epsg: u32,
178    datum: Datum,
179    method: ProjectionMethod,
180    linear_unit: LinearUnit,
181    name: &'static str,
182}
183
184impl ProjectedCrsDef {
185    pub const fn new(
186        epsg: u32,
187        datum: Datum,
188        method: ProjectionMethod,
189        linear_unit: LinearUnit,
190        name: &'static str,
191    ) -> Self {
192        Self::new_with_base_geographic_crs(epsg, 0, datum, method, linear_unit, name)
193    }
194
195    pub const fn new_with_base_geographic_crs(
196        epsg: u32,
197        base_geographic_crs_epsg: u32,
198        datum: Datum,
199        method: ProjectionMethod,
200        linear_unit: LinearUnit,
201        name: &'static str,
202    ) -> Self {
203        Self {
204            epsg,
205            base_geographic_crs_epsg,
206            datum,
207            method,
208            linear_unit,
209            name,
210        }
211    }
212
213    pub const fn epsg(&self) -> u32 {
214        self.epsg
215    }
216
217    pub const fn datum(&self) -> &Datum {
218        &self.datum
219    }
220
221    pub const fn base_geographic_crs_epsg(&self) -> u32 {
222        self.base_geographic_crs_epsg
223    }
224
225    pub const fn method(&self) -> ProjectionMethod {
226        self.method
227    }
228
229    pub const fn linear_unit(&self) -> LinearUnit {
230        self.linear_unit
231    }
232
233    pub const fn linear_unit_to_meter(&self) -> f64 {
234        self.linear_unit.meters_per_unit()
235    }
236
237    pub const fn name(&self) -> &'static str {
238        self.name
239    }
240}
241
242/// All supported projection methods with their parameters.
243///
244/// Angle parameters are stored in **degrees**. Conversion to radians happens
245/// at projection construction time (once), not per-transform.
246#[derive(Debug, Clone, Copy, PartialEq)]
247pub enum ProjectionMethod {
248    /// Web Mercator (EPSG:3857) — spherical Mercator on WGS84 semi-major axis.
249    WebMercator,
250
251    /// Transverse Mercator (includes UTM zones).
252    TransverseMercator {
253        /// Central meridian (degrees).
254        lon0: f64,
255        /// Latitude of origin (degrees).
256        lat0: f64,
257        /// Scale factor on central meridian.
258        k0: f64,
259        /// False easting (meters).
260        false_easting: f64,
261        /// False northing (meters).
262        false_northing: f64,
263    },
264
265    /// Polar Stereographic.
266    PolarStereographic {
267        /// Central meridian / straight vertical longitude (degrees).
268        lon0: f64,
269        /// Latitude of true scale (degrees). Determines the hemisphere.
270        lat_ts: f64,
271        /// Scale factor (used when lat_ts = ±90°, otherwise derived from lat_ts).
272        k0: f64,
273        /// False easting (meters).
274        false_easting: f64,
275        /// False northing (meters).
276        false_northing: f64,
277    },
278
279    /// Lambert Conformal Conic (1SP or 2SP).
280    LambertConformalConic {
281        /// Central meridian (degrees).
282        lon0: f64,
283        /// Latitude of origin (degrees).
284        lat0: f64,
285        /// First standard parallel (degrees).
286        lat1: f64,
287        /// Second standard parallel (degrees). Set equal to lat1 for 1SP variant.
288        lat2: f64,
289        /// False easting (meters).
290        false_easting: f64,
291        /// False northing (meters).
292        false_northing: f64,
293    },
294
295    /// Albers Equal Area Conic.
296    AlbersEqualArea {
297        /// Central meridian (degrees).
298        lon0: f64,
299        /// Latitude of origin (degrees).
300        lat0: f64,
301        /// First standard parallel (degrees).
302        lat1: f64,
303        /// Second standard parallel (degrees).
304        lat2: f64,
305        /// False easting (meters).
306        false_easting: f64,
307        /// False northing (meters).
308        false_northing: f64,
309    },
310
311    /// Lambert Azimuthal Equal Area.
312    LambertAzimuthalEqualArea {
313        /// Longitude of natural origin (degrees).
314        lon0: f64,
315        /// Latitude of natural origin (degrees).
316        lat0: f64,
317        /// False easting (meters).
318        false_easting: f64,
319        /// False northing (meters).
320        false_northing: f64,
321    },
322
323    /// Lambert Azimuthal Equal Area (spherical).
324    LambertAzimuthalEqualAreaSpherical {
325        /// Longitude of natural origin (degrees).
326        lon0: f64,
327        /// Latitude of natural origin (degrees).
328        lat0: f64,
329        /// False easting (meters).
330        false_easting: f64,
331        /// False northing (meters).
332        false_northing: f64,
333    },
334
335    /// EPSG Oblique Stereographic (Roussilhe / double stereographic).
336    ObliqueStereographic {
337        /// Longitude of natural origin (degrees).
338        lon0: f64,
339        /// Latitude of natural origin (degrees).
340        lat0: f64,
341        /// Scale factor at natural origin.
342        k0: f64,
343        /// False easting (meters).
344        false_easting: f64,
345        /// False northing (meters).
346        false_northing: f64,
347    },
348
349    /// Hotine Oblique Mercator / Rectified Skew Orthomorphic.
350    HotineObliqueMercator {
351        /// Latitude of projection centre (degrees).
352        latc: f64,
353        /// Longitude of projection centre (degrees).
354        lonc: f64,
355        /// Azimuth of central line at projection centre (degrees clockwise from north).
356        azimuth: f64,
357        /// Angle from rectified to skew grid (degrees).
358        rectified_grid_angle: f64,
359        /// Scale factor at projection centre.
360        k0: f64,
361        /// False easting or easting at projection centre (meters).
362        false_easting: f64,
363        /// False northing or northing at projection centre (meters).
364        false_northing: f64,
365        /// EPSG variant B offsets the natural origin to the projection centre.
366        variant_b: bool,
367    },
368
369    /// Cassini-Soldner.
370    CassiniSoldner {
371        /// Longitude of natural origin (degrees).
372        lon0: f64,
373        /// Latitude of natural origin (degrees).
374        lat0: f64,
375        /// False easting (meters).
376        false_easting: f64,
377        /// False northing (meters).
378        false_northing: f64,
379    },
380
381    /// Standard Mercator (ellipsoidal, distinct from Web Mercator).
382    Mercator {
383        /// Central meridian (degrees).
384        lon0: f64,
385        /// Latitude of true scale (degrees). 0 for 1SP variant.
386        lat_ts: f64,
387        /// Scale factor (for 1SP when lat_ts = 0).
388        k0: f64,
389        /// False easting (meters).
390        false_easting: f64,
391        /// False northing (meters).
392        false_northing: f64,
393    },
394
395    /// Equidistant Cylindrical / Plate Carrée.
396    EquidistantCylindrical {
397        /// Central meridian (degrees).
398        lon0: f64,
399        /// Latitude of true scale (degrees).
400        lat_ts: f64,
401        /// False easting (meters).
402        false_easting: f64,
403        /// False northing (meters).
404        false_northing: f64,
405    },
406}
407
408fn projection_methods_equivalent(a: &ProjectionMethod, b: &ProjectionMethod) -> bool {
409    match (a, b) {
410        (ProjectionMethod::WebMercator, ProjectionMethod::WebMercator) => true,
411        (
412            ProjectionMethod::TransverseMercator {
413                lon0: a_lon0,
414                lat0: a_lat0,
415                k0: a_k0,
416                false_easting: a_false_easting,
417                false_northing: a_false_northing,
418            },
419            ProjectionMethod::TransverseMercator {
420                lon0: b_lon0,
421                lat0: b_lat0,
422                k0: b_k0,
423                false_easting: b_false_easting,
424                false_northing: b_false_northing,
425            },
426        ) => {
427            approx_eq(*a_lon0, *b_lon0)
428                && approx_eq(*a_lat0, *b_lat0)
429                && approx_eq(*a_k0, *b_k0)
430                && approx_eq(*a_false_easting, *b_false_easting)
431                && approx_eq(*a_false_northing, *b_false_northing)
432        }
433        (
434            ProjectionMethod::PolarStereographic {
435                lon0: a_lon0,
436                lat_ts: a_lat_ts,
437                k0: a_k0,
438                false_easting: a_false_easting,
439                false_northing: a_false_northing,
440            },
441            ProjectionMethod::PolarStereographic {
442                lon0: b_lon0,
443                lat_ts: b_lat_ts,
444                k0: b_k0,
445                false_easting: b_false_easting,
446                false_northing: b_false_northing,
447            },
448        ) => {
449            approx_eq(*a_lon0, *b_lon0)
450                && approx_eq(*a_lat_ts, *b_lat_ts)
451                && approx_eq(*a_k0, *b_k0)
452                && approx_eq(*a_false_easting, *b_false_easting)
453                && approx_eq(*a_false_northing, *b_false_northing)
454        }
455        (
456            ProjectionMethod::LambertConformalConic {
457                lon0: a_lon0,
458                lat0: a_lat0,
459                lat1: a_lat1,
460                lat2: a_lat2,
461                false_easting: a_false_easting,
462                false_northing: a_false_northing,
463            },
464            ProjectionMethod::LambertConformalConic {
465                lon0: b_lon0,
466                lat0: b_lat0,
467                lat1: b_lat1,
468                lat2: b_lat2,
469                false_easting: b_false_easting,
470                false_northing: b_false_northing,
471            },
472        ) => {
473            approx_eq(*a_lon0, *b_lon0)
474                && approx_eq(*a_lat0, *b_lat0)
475                && approx_eq(*a_lat1, *b_lat1)
476                && approx_eq(*a_lat2, *b_lat2)
477                && approx_eq(*a_false_easting, *b_false_easting)
478                && approx_eq(*a_false_northing, *b_false_northing)
479        }
480        (
481            ProjectionMethod::AlbersEqualArea {
482                lon0: a_lon0,
483                lat0: a_lat0,
484                lat1: a_lat1,
485                lat2: a_lat2,
486                false_easting: a_false_easting,
487                false_northing: a_false_northing,
488            },
489            ProjectionMethod::AlbersEqualArea {
490                lon0: b_lon0,
491                lat0: b_lat0,
492                lat1: b_lat1,
493                lat2: b_lat2,
494                false_easting: b_false_easting,
495                false_northing: b_false_northing,
496            },
497        ) => {
498            approx_eq(*a_lon0, *b_lon0)
499                && approx_eq(*a_lat0, *b_lat0)
500                && approx_eq(*a_lat1, *b_lat1)
501                && approx_eq(*a_lat2, *b_lat2)
502                && approx_eq(*a_false_easting, *b_false_easting)
503                && approx_eq(*a_false_northing, *b_false_northing)
504        }
505        (
506            ProjectionMethod::LambertAzimuthalEqualArea {
507                lon0: a_lon0,
508                lat0: a_lat0,
509                false_easting: a_false_easting,
510                false_northing: a_false_northing,
511            },
512            ProjectionMethod::LambertAzimuthalEqualArea {
513                lon0: b_lon0,
514                lat0: b_lat0,
515                false_easting: b_false_easting,
516                false_northing: b_false_northing,
517            },
518        ) => {
519            approx_eq(*a_lon0, *b_lon0)
520                && approx_eq(*a_lat0, *b_lat0)
521                && approx_eq(*a_false_easting, *b_false_easting)
522                && approx_eq(*a_false_northing, *b_false_northing)
523        }
524        (
525            ProjectionMethod::LambertAzimuthalEqualAreaSpherical {
526                lon0: a_lon0,
527                lat0: a_lat0,
528                false_easting: a_false_easting,
529                false_northing: a_false_northing,
530            },
531            ProjectionMethod::LambertAzimuthalEqualAreaSpherical {
532                lon0: b_lon0,
533                lat0: b_lat0,
534                false_easting: b_false_easting,
535                false_northing: b_false_northing,
536            },
537        ) => {
538            approx_eq(*a_lon0, *b_lon0)
539                && approx_eq(*a_lat0, *b_lat0)
540                && approx_eq(*a_false_easting, *b_false_easting)
541                && approx_eq(*a_false_northing, *b_false_northing)
542        }
543        (
544            ProjectionMethod::ObliqueStereographic {
545                lon0: a_lon0,
546                lat0: a_lat0,
547                k0: a_k0,
548                false_easting: a_false_easting,
549                false_northing: a_false_northing,
550            },
551            ProjectionMethod::ObliqueStereographic {
552                lon0: b_lon0,
553                lat0: b_lat0,
554                k0: b_k0,
555                false_easting: b_false_easting,
556                false_northing: b_false_northing,
557            },
558        ) => {
559            approx_eq(*a_lon0, *b_lon0)
560                && approx_eq(*a_lat0, *b_lat0)
561                && approx_eq(*a_k0, *b_k0)
562                && approx_eq(*a_false_easting, *b_false_easting)
563                && approx_eq(*a_false_northing, *b_false_northing)
564        }
565        (
566            ProjectionMethod::HotineObliqueMercator {
567                latc: a_latc,
568                lonc: a_lonc,
569                azimuth: a_azimuth,
570                rectified_grid_angle: a_rectified_grid_angle,
571                k0: a_k0,
572                false_easting: a_false_easting,
573                false_northing: a_false_northing,
574                variant_b: a_variant_b,
575            },
576            ProjectionMethod::HotineObliqueMercator {
577                latc: b_latc,
578                lonc: b_lonc,
579                azimuth: b_azimuth,
580                rectified_grid_angle: b_rectified_grid_angle,
581                k0: b_k0,
582                false_easting: b_false_easting,
583                false_northing: b_false_northing,
584                variant_b: b_variant_b,
585            },
586        ) => {
587            a_variant_b == b_variant_b
588                && approx_eq(*a_latc, *b_latc)
589                && approx_eq(*a_lonc, *b_lonc)
590                && approx_eq(*a_azimuth, *b_azimuth)
591                && approx_eq(*a_rectified_grid_angle, *b_rectified_grid_angle)
592                && approx_eq(*a_k0, *b_k0)
593                && approx_eq(*a_false_easting, *b_false_easting)
594                && approx_eq(*a_false_northing, *b_false_northing)
595        }
596        (
597            ProjectionMethod::CassiniSoldner {
598                lon0: a_lon0,
599                lat0: a_lat0,
600                false_easting: a_false_easting,
601                false_northing: a_false_northing,
602            },
603            ProjectionMethod::CassiniSoldner {
604                lon0: b_lon0,
605                lat0: b_lat0,
606                false_easting: b_false_easting,
607                false_northing: b_false_northing,
608            },
609        ) => {
610            approx_eq(*a_lon0, *b_lon0)
611                && approx_eq(*a_lat0, *b_lat0)
612                && approx_eq(*a_false_easting, *b_false_easting)
613                && approx_eq(*a_false_northing, *b_false_northing)
614        }
615        (
616            ProjectionMethod::Mercator {
617                lon0: a_lon0,
618                lat_ts: a_lat_ts,
619                k0: a_k0,
620                false_easting: a_false_easting,
621                false_northing: a_false_northing,
622            },
623            ProjectionMethod::Mercator {
624                lon0: b_lon0,
625                lat_ts: b_lat_ts,
626                k0: b_k0,
627                false_easting: b_false_easting,
628                false_northing: b_false_northing,
629            },
630        ) => {
631            approx_eq(*a_lon0, *b_lon0)
632                && approx_eq(*a_lat_ts, *b_lat_ts)
633                && approx_eq(*a_k0, *b_k0)
634                && approx_eq(*a_false_easting, *b_false_easting)
635                && approx_eq(*a_false_northing, *b_false_northing)
636        }
637        (
638            ProjectionMethod::EquidistantCylindrical {
639                lon0: a_lon0,
640                lat_ts: a_lat_ts,
641                false_easting: a_false_easting,
642                false_northing: a_false_northing,
643            },
644            ProjectionMethod::EquidistantCylindrical {
645                lon0: b_lon0,
646                lat_ts: b_lat_ts,
647                false_easting: b_false_easting,
648                false_northing: b_false_northing,
649            },
650        ) => {
651            approx_eq(*a_lon0, *b_lon0)
652                && approx_eq(*a_lat_ts, *b_lat_ts)
653                && approx_eq(*a_false_easting, *b_false_easting)
654                && approx_eq(*a_false_northing, *b_false_northing)
655        }
656        _ => false,
657    }
658}
659
660fn approx_eq(a: f64, b: f64) -> bool {
661    (a - b).abs() < 1e-12
662}
663
664#[cfg(test)]
665mod tests {
666    use super::*;
667    use crate::datum;
668
669    #[test]
670    fn geographic_crs_is_geographic() {
671        let crs = CrsDef::Geographic(GeographicCrsDef::new(4326, datum::WGS84, "WGS 84"));
672        assert!(crs.is_geographic());
673        assert!(!crs.is_projected());
674        assert_eq!(crs.epsg(), 4326);
675    }
676
677    #[test]
678    fn projected_crs_is_projected() {
679        let crs = CrsDef::Projected(ProjectedCrsDef::new(
680            3857,
681            datum::WGS84,
682            ProjectionMethod::WebMercator,
683            LinearUnit::metre(),
684            "WGS 84 / Pseudo-Mercator",
685        ));
686        assert!(crs.is_projected());
687        assert!(!crs.is_geographic());
688        assert_eq!(crs.epsg(), 3857);
689    }
690
691    #[test]
692    fn linear_unit_validates_positive_finite_conversion() {
693        assert!(LinearUnit::from_meters_per_unit(0.3048).is_ok());
694        assert!(LinearUnit::from_meters_per_unit(0.0).is_err());
695        assert!(LinearUnit::from_meters_per_unit(f64::NAN).is_err());
696    }
697}