Skip to main content

proj_core/
transform.rs

1use crate::coord::{Bounds, Coord, Coord3D, Transformable, Transformable3D};
2use crate::crs::{CrsDef, ProjectionMethod};
3use crate::datum::Datum;
4use crate::error::{Error, Result};
5use crate::geocentric;
6use crate::helmert;
7use crate::projection::{make_projection, ProjectionImpl};
8use crate::registry;
9
10/// A reusable coordinate transformation between two CRS.
11///
12/// Create once with [`Transform::new`], then call [`convert`](Transform::convert)
13/// or [`convert_3d`](Transform::convert_3d) for each coordinate. All input/output
14/// coordinates use the CRS's native units: degrees (lon/lat) for geographic CRS,
15/// meters for projected CRS. For 3D coordinates, the third ordinate is preserved
16/// unchanged by the current horizontal-only pipeline.
17///
18/// # Example
19///
20/// ```
21/// use proj_core::Transform;
22///
23/// let t = Transform::new("EPSG:4326", "EPSG:3857").unwrap();
24/// let (x, y) = t.convert((-74.006, 40.7128)).unwrap();
25/// assert!((x - (-8238310.0)).abs() < 1.0);
26/// ```
27pub struct Transform {
28    source: CrsDef,
29    target: CrsDef,
30    pipeline: TransformPipeline,
31}
32
33enum TransformPipeline {
34    /// Source and target are the same CRS — no-op.
35    Identity,
36    /// Same datum, geographic source → projected target.
37    SameDatumForward { forward: Box<dyn ProjectionImpl> },
38    /// Same datum, projected source → geographic target.
39    SameDatumInverse { inverse: Box<dyn ProjectionImpl> },
40    /// Same datum, projected source → projected target.
41    SameDatumBoth {
42        inverse: Box<dyn ProjectionImpl>,
43        forward: Box<dyn ProjectionImpl>,
44    },
45    /// Different datums — full pipeline with geocentric conversion and Helmert shift.
46    DatumShift {
47        inverse: Option<Box<dyn ProjectionImpl>>,
48        source_datum: Datum,
49        target_datum: Datum,
50        forward: Option<Box<dyn ProjectionImpl>>,
51    },
52}
53
54impl Transform {
55    /// Create a transform from authority code strings (e.g., `"EPSG:4326"`).
56    ///
57    /// This is the primary constructor, matching the API pattern from C PROJ's
58    /// `Proj::new_known_crs()`.
59    pub fn new(from_crs: &str, to_crs: &str) -> Result<Self> {
60        let source = registry::lookup_authority_code(from_crs)?;
61        let target = registry::lookup_authority_code(to_crs)?;
62        Self::from_crs_defs(&source, &target)
63    }
64
65    /// Create a transform from EPSG codes directly.
66    pub fn from_epsg(from: u32, to: u32) -> Result<Self> {
67        let source = registry::lookup_epsg(from)
68            .ok_or_else(|| Error::UnknownCrs(format!("unknown EPSG code: {from}")))?;
69        let target = registry::lookup_epsg(to)
70            .ok_or_else(|| Error::UnknownCrs(format!("unknown EPSG code: {to}")))?;
71        Self::from_crs_defs(&source, &target)
72    }
73
74    /// Create a transform from explicit CRS definitions.
75    ///
76    /// Use this for custom CRS not in the built-in registry.
77    pub fn from_crs_defs(from: &CrsDef, to: &CrsDef) -> Result<Self> {
78        // Identity check for known EPSG codes and semantically identical custom CRS definitions.
79        if (from.epsg() != 0 && from.epsg() == to.epsg()) || same_crs_definition(from, to) {
80            return Ok(Self {
81                source: *from,
82                target: *to,
83                pipeline: TransformPipeline::Identity,
84            });
85        }
86
87        let source_datum = from.datum();
88        let target_datum = to.datum();
89
90        let same_datum = source_datum.same_datum(target_datum)
91            || (source_datum.is_wgs84_compatible() && target_datum.is_wgs84_compatible());
92
93        let pipeline = if same_datum {
94            // Optimized path: no datum shift needed
95            match (from, to) {
96                (CrsDef::Geographic(_), CrsDef::Geographic(_)) => TransformPipeline::Identity,
97                (CrsDef::Geographic(_), CrsDef::Projected(p)) => {
98                    let forward = make_projection(&p.method, &p.datum)?;
99                    TransformPipeline::SameDatumForward { forward }
100                }
101                (CrsDef::Projected(p), CrsDef::Geographic(_)) => {
102                    let inverse = make_projection(&p.method, &p.datum)?;
103                    TransformPipeline::SameDatumInverse { inverse }
104                }
105                (CrsDef::Projected(p_from), CrsDef::Projected(p_to)) => {
106                    let inverse = make_projection(&p_from.method, &p_from.datum)?;
107                    let forward = make_projection(&p_to.method, &p_to.datum)?;
108                    TransformPipeline::SameDatumBoth { inverse, forward }
109                }
110            }
111        } else {
112            // Both datums must have a path to WGS84
113            if !source_datum.is_wgs84_compatible() && source_datum.to_wgs84.is_none() {
114                return Err(Error::UnsupportedProjection(format!(
115                    "source CRS EPSG:{} has no known datum shift to WGS84",
116                    from.epsg()
117                )));
118            }
119            if !target_datum.is_wgs84_compatible() && target_datum.to_wgs84.is_none() {
120                return Err(Error::UnsupportedProjection(format!(
121                    "target CRS EPSG:{} has no known datum shift to WGS84",
122                    to.epsg()
123                )));
124            }
125
126            let inverse = match from {
127                CrsDef::Projected(p) => Some(make_projection(&p.method, &p.datum)?),
128                CrsDef::Geographic(_) => None,
129            };
130            let forward = match to {
131                CrsDef::Projected(p) => Some(make_projection(&p.method, &p.datum)?),
132                CrsDef::Geographic(_) => None,
133            };
134
135            TransformPipeline::DatumShift {
136                inverse,
137                source_datum: *source_datum,
138                target_datum: *target_datum,
139                forward,
140            }
141        };
142
143        Ok(Self {
144            source: *from,
145            target: *to,
146            pipeline,
147        })
148    }
149
150    /// Transform a single coordinate.
151    ///
152    /// Input and output units are the native units of the respective CRS:
153    /// degrees for geographic CRS, meters for projected CRS.
154    ///
155    /// The return type matches the input type:
156    /// - `(f64, f64)` in → `(f64, f64)` out
157    /// - `Coord` in → `Coord` out
158    /// - `geo_types::Coord<f64>` in → `geo_types::Coord<f64>` out (with `geo-types` feature)
159    pub fn convert<T: Transformable>(&self, coord: T) -> Result<T> {
160        let c = coord.into_coord();
161        let result = self.convert_coord(c)?;
162        Ok(T::from_coord(result))
163    }
164
165    /// Transform a single 3D coordinate.
166    ///
167    /// The horizontal components use the CRS's native units:
168    /// degrees for geographic CRS and meters for projected CRS.
169    /// The vertical component is carried through unchanged because the current
170    /// CRS model is horizontal-only.
171    ///
172    /// The return type matches the input type:
173    /// - `(f64, f64, f64)` in → `(f64, f64, f64)` out
174    /// - `Coord3D` in → `Coord3D` out
175    pub fn convert_3d<T: Transformable3D>(&self, coord: T) -> Result<T> {
176        let c = coord.into_coord3d();
177        let result = self.convert_coord3d(c)?;
178        Ok(T::from_coord3d(result))
179    }
180
181    /// Return the source CRS definition for this transform.
182    pub fn source_crs(&self) -> &CrsDef {
183        &self.source
184    }
185
186    /// Return the target CRS definition for this transform.
187    pub fn target_crs(&self) -> &CrsDef {
188        &self.target
189    }
190
191    /// Build the inverse transform by swapping the source and target CRS.
192    pub fn inverse(&self) -> Result<Self> {
193        Self::from_crs_defs(&self.target, &self.source)
194    }
195
196    /// Reproject a 2D bounding box by sampling its perimeter.
197    ///
198    /// `densify_points` controls how many evenly spaced interior points are sampled
199    /// on each edge between the corners. `0` samples only the four corners.
200    ///
201    /// The returned bounds are axis-aligned in the target CRS. Geographic outputs
202    /// that cross the antimeridian are not normalized into a wrapped representation.
203    pub fn transform_bounds(&self, bounds: Bounds, densify_points: usize) -> Result<Bounds> {
204        if !bounds.is_valid() {
205            return Err(Error::OutOfRange(
206                "bounds must be finite and satisfy min <= max".into(),
207            ));
208        }
209
210        let segments = densify_points
211            .checked_add(1)
212            .ok_or_else(|| Error::OutOfRange("densify point count is too large".into()))?;
213
214        let mut transformed: Option<Bounds> = None;
215        for i in 0..=segments {
216            let t = i as f64 / segments as f64;
217            let x = bounds.min_x + bounds.width() * t;
218            let y = bounds.min_y + bounds.height() * t;
219
220            for sample in [
221                Coord::new(x, bounds.min_y),
222                Coord::new(x, bounds.max_y),
223                Coord::new(bounds.min_x, y),
224                Coord::new(bounds.max_x, y),
225            ] {
226                let coord = self.convert_coord(sample)?;
227                if let Some(accum) = &mut transformed {
228                    accum.expand_to_include(coord);
229                } else {
230                    transformed = Some(Bounds::new(coord.x, coord.y, coord.x, coord.y));
231                }
232            }
233        }
234
235        transformed.ok_or_else(|| Error::OutOfRange("failed to sample bounds".into()))
236    }
237
238    /// Transform a single `Coord` value.
239    fn convert_coord(&self, c: Coord) -> Result<Coord> {
240        let result = self.convert_coord3d(Coord3D::new(c.x, c.y, 0.0))?;
241        Ok(Coord {
242            x: result.x,
243            y: result.y,
244        })
245    }
246
247    /// Transform a single `Coord3D` value.
248    fn convert_coord3d(&self, c: Coord3D) -> Result<Coord3D> {
249        match &self.pipeline {
250            TransformPipeline::Identity => Ok(c),
251
252            TransformPipeline::SameDatumForward { forward } => {
253                // Source is geographic (degrees) → radians → forward project → meters
254                let lon_rad = c.x.to_radians();
255                let lat_rad = c.y.to_radians();
256                let (x, y) = forward.forward(lon_rad, lat_rad)?;
257                Ok(Coord3D { x, y, z: c.z })
258            }
259
260            TransformPipeline::SameDatumInverse { inverse } => {
261                // Source is projected (meters) → inverse project → radians → degrees
262                let (lon_rad, lat_rad) = inverse.inverse(c.x, c.y)?;
263                Ok(Coord3D {
264                    x: lon_rad.to_degrees(),
265                    y: lat_rad.to_degrees(),
266                    z: c.z,
267                })
268            }
269
270            TransformPipeline::SameDatumBoth { inverse, forward } => {
271                // Projected → inverse → radians → forward → projected
272                let (lon_rad, lat_rad) = inverse.inverse(c.x, c.y)?;
273                let (x, y) = forward.forward(lon_rad, lat_rad)?;
274                Ok(Coord3D { x, y, z: c.z })
275            }
276
277            TransformPipeline::DatumShift {
278                inverse,
279                source_datum,
280                target_datum,
281                forward,
282            } => {
283                // Step 1: Get geographic coords in radians on source datum
284                let (lon_rad, lat_rad) = if let Some(inv) = inverse {
285                    inv.inverse(c.x, c.y)?
286                } else {
287                    (c.x.to_radians(), c.y.to_radians())
288                };
289
290                // Step 2: Source geodetic → geocentric (source ellipsoid)
291                let (x, y, z) = geocentric::geodetic_to_geocentric(
292                    &source_datum.ellipsoid,
293                    lon_rad,
294                    lat_rad,
295                    0.0,
296                );
297
298                // Step 3: Helmert: source datum → WGS84
299                let (x2, y2, z2) = if let Some(params) = &source_datum.to_wgs84 {
300                    helmert::helmert_forward(params, x, y, z)
301                } else {
302                    (x, y, z) // already WGS84
303                };
304
305                // Step 4: Helmert: WGS84 → target datum (inverse)
306                let (x3, y3, z3) = if let Some(params) = &target_datum.to_wgs84 {
307                    helmert::helmert_inverse(params, x2, y2, z2)
308                } else {
309                    (x2, y2, z2) // target is WGS84
310                };
311
312                // Step 5: Geocentric → geodetic (target ellipsoid)
313                let (lon_out, lat_out, _h_out) =
314                    geocentric::geocentric_to_geodetic(&target_datum.ellipsoid, x3, y3, z3);
315
316                // Step 6: Forward project if target is projected, else convert to degrees
317                if let Some(fwd) = forward {
318                    let (x, y) = fwd.forward(lon_out, lat_out)?;
319                    Ok(Coord3D { x, y, z: c.z })
320                } else {
321                    Ok(Coord3D {
322                        x: lon_out.to_degrees(),
323                        y: lat_out.to_degrees(),
324                        z: c.z,
325                    })
326                }
327            }
328        }
329    }
330
331    /// Batch transform (sequential).
332    pub fn convert_batch<T: Transformable + Clone>(&self, coords: &[T]) -> Result<Vec<T>> {
333        coords.iter().map(|c| self.convert(c.clone())).collect()
334    }
335
336    /// Batch transform of 3D coordinates (sequential).
337    pub fn convert_batch_3d<T: Transformable3D + Clone>(&self, coords: &[T]) -> Result<Vec<T>> {
338        coords.iter().map(|c| self.convert_3d(c.clone())).collect()
339    }
340
341    /// Batch transform with Rayon parallelism.
342    #[cfg(feature = "rayon")]
343    pub fn convert_batch_parallel<T: Transformable + Send + Sync + Clone>(
344        &self,
345        coords: &[T],
346    ) -> Result<Vec<T>> {
347        use rayon::prelude::*;
348        coords.par_iter().map(|c| self.convert(c.clone())).collect()
349    }
350
351    /// Batch transform of 3D coordinates with Rayon parallelism.
352    #[cfg(feature = "rayon")]
353    pub fn convert_batch_parallel_3d<T: Transformable3D + Send + Sync + Clone>(
354        &self,
355        coords: &[T],
356    ) -> Result<Vec<T>> {
357        use rayon::prelude::*;
358        coords
359            .par_iter()
360            .map(|c| self.convert_3d(c.clone()))
361            .collect()
362    }
363}
364
365fn same_crs_definition(from: &CrsDef, to: &CrsDef) -> bool {
366    match (from, to) {
367        (CrsDef::Geographic(a), CrsDef::Geographic(b)) => a.datum.same_datum(&b.datum),
368        (CrsDef::Projected(a), CrsDef::Projected(b)) => {
369            a.datum.same_datum(&b.datum) && projection_methods_equivalent(&a.method, &b.method)
370        }
371        _ => false,
372    }
373}
374
375fn projection_methods_equivalent(a: &ProjectionMethod, b: &ProjectionMethod) -> bool {
376    match (a, b) {
377        (ProjectionMethod::WebMercator, ProjectionMethod::WebMercator) => true,
378        (
379            ProjectionMethod::TransverseMercator {
380                lon0: a_lon0,
381                lat0: a_lat0,
382                k0: a_k0,
383                false_easting: a_false_easting,
384                false_northing: a_false_northing,
385            },
386            ProjectionMethod::TransverseMercator {
387                lon0: b_lon0,
388                lat0: b_lat0,
389                k0: b_k0,
390                false_easting: b_false_easting,
391                false_northing: b_false_northing,
392            },
393        ) => {
394            approx_eq(*a_lon0, *b_lon0)
395                && approx_eq(*a_lat0, *b_lat0)
396                && approx_eq(*a_k0, *b_k0)
397                && approx_eq(*a_false_easting, *b_false_easting)
398                && approx_eq(*a_false_northing, *b_false_northing)
399        }
400        (
401            ProjectionMethod::PolarStereographic {
402                lon0: a_lon0,
403                lat_ts: a_lat_ts,
404                k0: a_k0,
405                false_easting: a_false_easting,
406                false_northing: a_false_northing,
407            },
408            ProjectionMethod::PolarStereographic {
409                lon0: b_lon0,
410                lat_ts: b_lat_ts,
411                k0: b_k0,
412                false_easting: b_false_easting,
413                false_northing: b_false_northing,
414            },
415        ) => {
416            approx_eq(*a_lon0, *b_lon0)
417                && approx_eq(*a_lat_ts, *b_lat_ts)
418                && approx_eq(*a_k0, *b_k0)
419                && approx_eq(*a_false_easting, *b_false_easting)
420                && approx_eq(*a_false_northing, *b_false_northing)
421        }
422        (
423            ProjectionMethod::LambertConformalConic {
424                lon0: a_lon0,
425                lat0: a_lat0,
426                lat1: a_lat1,
427                lat2: a_lat2,
428                false_easting: a_false_easting,
429                false_northing: a_false_northing,
430            },
431            ProjectionMethod::LambertConformalConic {
432                lon0: b_lon0,
433                lat0: b_lat0,
434                lat1: b_lat1,
435                lat2: b_lat2,
436                false_easting: b_false_easting,
437                false_northing: b_false_northing,
438            },
439        ) => {
440            approx_eq(*a_lon0, *b_lon0)
441                && approx_eq(*a_lat0, *b_lat0)
442                && approx_eq(*a_lat1, *b_lat1)
443                && approx_eq(*a_lat2, *b_lat2)
444                && approx_eq(*a_false_easting, *b_false_easting)
445                && approx_eq(*a_false_northing, *b_false_northing)
446        }
447        (
448            ProjectionMethod::AlbersEqualArea {
449                lon0: a_lon0,
450                lat0: a_lat0,
451                lat1: a_lat1,
452                lat2: a_lat2,
453                false_easting: a_false_easting,
454                false_northing: a_false_northing,
455            },
456            ProjectionMethod::AlbersEqualArea {
457                lon0: b_lon0,
458                lat0: b_lat0,
459                lat1: b_lat1,
460                lat2: b_lat2,
461                false_easting: b_false_easting,
462                false_northing: b_false_northing,
463            },
464        ) => {
465            approx_eq(*a_lon0, *b_lon0)
466                && approx_eq(*a_lat0, *b_lat0)
467                && approx_eq(*a_lat1, *b_lat1)
468                && approx_eq(*a_lat2, *b_lat2)
469                && approx_eq(*a_false_easting, *b_false_easting)
470                && approx_eq(*a_false_northing, *b_false_northing)
471        }
472        (
473            ProjectionMethod::Mercator {
474                lon0: a_lon0,
475                lat_ts: a_lat_ts,
476                k0: a_k0,
477                false_easting: a_false_easting,
478                false_northing: a_false_northing,
479            },
480            ProjectionMethod::Mercator {
481                lon0: b_lon0,
482                lat_ts: b_lat_ts,
483                k0: b_k0,
484                false_easting: b_false_easting,
485                false_northing: b_false_northing,
486            },
487        ) => {
488            approx_eq(*a_lon0, *b_lon0)
489                && approx_eq(*a_lat_ts, *b_lat_ts)
490                && approx_eq(*a_k0, *b_k0)
491                && approx_eq(*a_false_easting, *b_false_easting)
492                && approx_eq(*a_false_northing, *b_false_northing)
493        }
494        (
495            ProjectionMethod::EquidistantCylindrical {
496                lon0: a_lon0,
497                lat_ts: a_lat_ts,
498                false_easting: a_false_easting,
499                false_northing: a_false_northing,
500            },
501            ProjectionMethod::EquidistantCylindrical {
502                lon0: b_lon0,
503                lat_ts: b_lat_ts,
504                false_easting: b_false_easting,
505                false_northing: b_false_northing,
506            },
507        ) => {
508            approx_eq(*a_lon0, *b_lon0)
509                && approx_eq(*a_lat_ts, *b_lat_ts)
510                && approx_eq(*a_false_easting, *b_false_easting)
511                && approx_eq(*a_false_northing, *b_false_northing)
512        }
513        _ => false,
514    }
515}
516
517fn approx_eq(a: f64, b: f64) -> bool {
518    (a - b).abs() < 1e-12
519}
520
521#[cfg(test)]
522mod tests {
523    use super::*;
524    use crate::crs::{CrsDef, ProjectedCrsDef, ProjectionMethod};
525    use crate::datum;
526
527    #[test]
528    fn identity_same_crs() {
529        let t = Transform::new("EPSG:4326", "EPSG:4326").unwrap();
530        let (x, y) = t.convert((-74.006, 40.7128)).unwrap();
531        assert_eq!(x, -74.006);
532        assert_eq!(y, 40.7128);
533    }
534
535    #[test]
536    fn wgs84_to_web_mercator() {
537        let t = Transform::new("EPSG:4326", "EPSG:3857").unwrap();
538        let (x, y) = t.convert((-74.006, 40.7128)).unwrap();
539        // NYC in Web Mercator
540        assert!((x - (-8238310.0)).abs() < 100.0, "x = {x}");
541        assert!((y - 4970072.0).abs() < 100.0, "y = {y}");
542    }
543
544    #[test]
545    fn web_mercator_to_wgs84() {
546        let t = Transform::new("EPSG:3857", "EPSG:4326").unwrap();
547        let (lon, lat) = t.convert((-8238310.0, 4970072.0)).unwrap();
548        assert!((lon - (-74.006)).abs() < 0.001, "lon = {lon}");
549        assert!((lat - 40.7128).abs() < 0.001, "lat = {lat}");
550    }
551
552    #[test]
553    fn roundtrip_4326_3857() {
554        let fwd = Transform::new("EPSG:4326", "EPSG:3857").unwrap();
555        let inv = fwd.inverse().unwrap();
556
557        let original = (-74.0445, 40.6892);
558        let projected = fwd.convert(original).unwrap();
559        let back = inv.convert(projected).unwrap();
560
561        assert!(
562            (back.0 - original.0).abs() < 1e-8,
563            "lon: {} vs {}",
564            back.0,
565            original.0
566        );
567        assert!(
568            (back.1 - original.1).abs() < 1e-8,
569            "lat: {} vs {}",
570            back.1,
571            original.1
572        );
573    }
574
575    #[test]
576    fn wgs84_to_utm_18n() {
577        let t = Transform::new("EPSG:4326", "EPSG:32618").unwrap();
578        let (x, y) = t.convert((-74.006, 40.7128)).unwrap();
579        assert!((x - 583960.0).abs() < 1.0, "easting = {x}");
580        assert!(y > 4_500_000.0 && y < 4_510_000.0, "northing = {y}");
581    }
582
583    #[test]
584    fn utm_to_web_mercator() {
585        // Projected → Projected (SameDatumBoth)
586        let t = Transform::new("EPSG:32618", "EPSG:3857").unwrap();
587        // NYC in UTM 18N
588        let (x, _y) = t.convert((583960.0, 4507523.0)).unwrap();
589        // Should be near NYC in Web Mercator
590        assert!((x - (-8238310.0)).abs() < 200.0, "x = {x}");
591    }
592
593    #[test]
594    fn wgs84_to_polar_stereo_3413() {
595        let t = Transform::new("EPSG:4326", "EPSG:3413").unwrap();
596        let (x, y) = t.convert((-45.0, 90.0)).unwrap();
597        // North pole should be at origin for EPSG:3413
598        assert!(x.abs() < 1.0, "x = {x}");
599        assert!(y.abs() < 1.0, "y = {y}");
600    }
601
602    #[test]
603    fn roundtrip_4326_3413() {
604        let fwd = Transform::new("EPSG:4326", "EPSG:3413").unwrap();
605        let inv = fwd.inverse().unwrap();
606
607        let original = (-45.0, 75.0);
608        let projected = fwd.convert(original).unwrap();
609        let back = inv.convert(projected).unwrap();
610
611        assert!(
612            (back.0 - original.0).abs() < 1e-6,
613            "lon: {} vs {}",
614            back.0,
615            original.0
616        );
617        assert!(
618            (back.1 - original.1).abs() < 1e-6,
619            "lat: {} vs {}",
620            back.1,
621            original.1
622        );
623    }
624
625    #[test]
626    fn geographic_to_geographic_same_datum_is_identity() {
627        // NAD83 and WGS84 are both WGS84-compatible
628        let t = Transform::new("EPSG:4269", "EPSG:4326").unwrap();
629        let (lon, lat) = t.convert((-74.006, 40.7128)).unwrap();
630        assert_eq!(lon, -74.006);
631        assert_eq!(lat, 40.7128);
632    }
633
634    #[test]
635    fn unknown_crs_error() {
636        let result = Transform::new("EPSG:99999", "EPSG:4326");
637        assert!(result.is_err());
638    }
639
640    #[test]
641    fn cross_datum_nad27_to_wgs84() {
642        let t = Transform::new("EPSG:4267", "EPSG:4326").unwrap();
643        let (lon, lat) = t.convert((-90.0, 45.0)).unwrap();
644        // NAD27→WGS84 shift is small (tens of meters → fraction of degree)
645        assert!((lon - (-90.0)).abs() < 0.01, "lon = {lon}");
646        assert!((lat - 45.0).abs() < 0.01, "lat = {lat}");
647    }
648
649    #[test]
650    fn cross_datum_roundtrip_nad27() {
651        let fwd = Transform::new("EPSG:4267", "EPSG:4326").unwrap();
652        let inv = fwd.inverse().unwrap();
653        let original = (-90.0, 45.0);
654        let shifted = fwd.convert(original).unwrap();
655        let back = inv.convert(shifted).unwrap();
656        assert!(
657            (back.0 - original.0).abs() < 1e-6,
658            "lon: {} vs {}",
659            back.0,
660            original.0
661        );
662        assert!(
663            (back.1 - original.1).abs() < 1e-6,
664            "lat: {} vs {}",
665            back.1,
666            original.1
667        );
668    }
669
670    #[test]
671    fn cross_datum_osgb36_to_wgs84() {
672        let t = Transform::new("EPSG:4277", "EPSG:4326").unwrap();
673        let (lon, lat) = t.convert((-0.1278, 51.5074)).unwrap(); // London
674                                                                 // OSGB36→WGS84 shift is larger due to 7-parameter Helmert
675        assert!((lon - (-0.1278)).abs() < 0.01, "lon = {lon}");
676        assert!((lat - 51.5074).abs() < 0.01, "lat = {lat}");
677    }
678
679    #[test]
680    fn wgs84_to_web_mercator_3d_preserves_height() {
681        let t = Transform::new("EPSG:4326", "EPSG:3857").unwrap();
682        let (x, y, z) = t.convert_3d((-74.006, 40.7128, 123.45)).unwrap();
683        assert!((x - (-8238310.0)).abs() < 100.0, "x = {x}");
684        assert!((y - 4970072.0).abs() < 100.0, "y = {y}");
685        assert!((z - 123.45).abs() < 1e-12, "z = {z}");
686    }
687
688    #[test]
689    fn cross_datum_roundtrip_nad27_3d() {
690        let fwd = Transform::new("EPSG:4267", "EPSG:4326").unwrap();
691        let inv = fwd.inverse().unwrap();
692        let original = (-90.0, 45.0, 250.0);
693        let shifted = fwd.convert_3d(original).unwrap();
694        let back = inv.convert_3d(shifted).unwrap();
695        assert!(
696            (back.0 - original.0).abs() < 1e-6,
697            "lon: {} vs {}",
698            back.0,
699            original.0
700        );
701        assert!(
702            (back.1 - original.1).abs() < 1e-6,
703            "lat: {} vs {}",
704            back.1,
705            original.1
706        );
707        assert!(
708            (back.2 - original.2).abs() < 1e-12,
709            "h: {} vs {}",
710            back.2,
711            original.2
712        );
713    }
714
715    #[test]
716    fn identical_custom_projected_crs_is_identity() {
717        let from = CrsDef::Projected(ProjectedCrsDef {
718            epsg: 0,
719            datum: datum::WGS84,
720            method: ProjectionMethod::WebMercator,
721            name: "Custom Web Mercator A",
722        });
723        let to = CrsDef::Projected(ProjectedCrsDef {
724            epsg: 0,
725            datum: datum::WGS84,
726            method: ProjectionMethod::WebMercator,
727            name: "Custom Web Mercator B",
728        });
729
730        let t = Transform::from_crs_defs(&from, &to).unwrap();
731        assert!(matches!(t.pipeline, TransformPipeline::Identity));
732    }
733
734    #[test]
735    fn inverse_exposes_swapped_crs() {
736        let fwd = Transform::new("EPSG:4326", "EPSG:3857").unwrap();
737        let inv = fwd.inverse().unwrap();
738
739        assert_eq!(fwd.source_crs().epsg(), 4326);
740        assert_eq!(fwd.target_crs().epsg(), 3857);
741        assert_eq!(inv.source_crs().epsg(), 3857);
742        assert_eq!(inv.target_crs().epsg(), 4326);
743    }
744
745    #[test]
746    fn transform_bounds_web_mercator() {
747        let t = Transform::new("EPSG:4326", "EPSG:3857").unwrap();
748        let bounds = Bounds::new(-74.3, 40.45, -73.65, 40.95);
749
750        let result = t.transform_bounds(bounds, 8).unwrap();
751
752        assert!(result.min_x < -8_200_000.0);
753        assert!(result.max_x < -8_100_000.0);
754        assert!(result.min_y > 4_900_000.0);
755        assert!(result.max_y > result.min_y);
756    }
757
758    #[test]
759    fn transform_bounds_rejects_invalid_input() {
760        let t = Transform::new("EPSG:4326", "EPSG:3857").unwrap();
761        let err = t
762            .transform_bounds(Bounds::new(10.0, 5.0, -10.0, 20.0), 0)
763            .unwrap_err();
764
765        assert!(matches!(err, Error::OutOfRange(_)));
766    }
767
768    #[test]
769    fn batch_transform() {
770        let t = Transform::new("EPSG:4326", "EPSG:3857").unwrap();
771        let coords: Vec<(f64, f64)> = (0..10)
772            .map(|i| (-74.0 + i as f64 * 0.1, 40.0 + i as f64 * 0.1))
773            .collect();
774
775        let results = t.convert_batch(&coords).unwrap();
776        assert_eq!(results.len(), 10);
777        for (x, _y) in &results {
778            assert!(*x < 0.0); // all points are west of prime meridian
779        }
780    }
781
782    #[test]
783    fn batch_transform_3d() {
784        let t = Transform::new("EPSG:4326", "EPSG:3857").unwrap();
785        let coords: Vec<(f64, f64, f64)> = (0..10)
786            .map(|i| (-74.0 + i as f64 * 0.1, 40.0 + i as f64 * 0.1, i as f64))
787            .collect();
788
789        let results = t.convert_batch_3d(&coords).unwrap();
790        assert_eq!(results.len(), 10);
791        for (index, (x, _y, z)) in results.iter().enumerate() {
792            assert!(*x < 0.0);
793            assert!((*z - index as f64).abs() < 1e-12);
794        }
795    }
796
797    #[test]
798    fn coord_type() {
799        let t = Transform::new("EPSG:4326", "EPSG:3857").unwrap();
800        let c = Coord::new(-74.006, 40.7128);
801        let result = t.convert(c).unwrap();
802        assert!((result.x - (-8238310.0)).abs() < 100.0);
803    }
804
805    #[cfg(feature = "geo-types")]
806    #[test]
807    fn geo_types_coord() {
808        let t = Transform::new("EPSG:4326", "EPSG:3857").unwrap();
809        let c = geo_types::Coord {
810            x: -74.006,
811            y: 40.7128,
812        };
813        let result: geo_types::Coord<f64> = t.convert(c).unwrap();
814        assert!((result.x - (-8238310.0)).abs() < 100.0);
815    }
816
817    #[cfg(feature = "rayon")]
818    #[test]
819    fn parallel_batch_transform() {
820        let t = Transform::new("EPSG:4326", "EPSG:3857").unwrap();
821        let coords: Vec<(f64, f64)> = (0..100)
822            .map(|i| (-74.0 + i as f64 * 0.01, 40.0 + i as f64 * 0.01))
823            .collect();
824
825        let results = t.convert_batch_parallel(&coords).unwrap();
826        assert_eq!(results.len(), 100);
827    }
828
829    #[cfg(feature = "rayon")]
830    #[test]
831    fn parallel_batch_transform_3d() {
832        let t = Transform::new("EPSG:4326", "EPSG:3857").unwrap();
833        let coords: Vec<(f64, f64, f64)> = (0..100)
834            .map(|i| (-74.0 + i as f64 * 0.01, 40.0 + i as f64 * 0.01, i as f64))
835            .collect();
836
837        let results = t.convert_batch_parallel_3d(&coords).unwrap();
838        assert_eq!(results.len(), 100);
839    }
840}