1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
//! Exact distortion factors computed via automatic differentiation.
//!
//! [`FactorsExact`] is the AD-aware companion to the numeric [`crate::Factors`]
//! struct. While [`crate::Factors`] uses finite differences, [`FactorsExact`]
//! computes the exact Jacobian via forward-mode AD (`Dual1<2>`), giving
//! machine-precision scale factors and angular distortion.
/// Exact distortion factors computed via automatic differentiation.
///
/// Companion to the numeric [`crate::Factors`] struct (kept for PROJ parity).
/// Constructed via [`FactorsExact::from_jacobian`] after evaluating a
/// `ProjectGeneric` implementation with `Dual1<2>` inputs.
#[derive(Debug, Clone, Copy)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct FactorsExact {
/// Meridional scale factor `h = ‖∂(x,y)/∂φ‖ / a`.
pub meridian_scale: f64,
/// Parallel scale factor `k = ‖∂(x,y)/∂λ‖ / (a·cos φ)`.
pub parallel_scale: f64,
/// Areal scale factor `s = |det J| / (a²·cos φ)`.
pub areal_scale: f64,
/// Maximum angular distortion `ω = 2·arcsin(|h−k|/(h+k))` in radians.
pub angular_distortion: f64,
/// Meridian convergence γ in radians: angle from grid north to true north.
pub meridian_convergence: f64,
/// Partial derivative ∂x/∂λ (longitude).
pub dx_dlam: f64,
/// Partial derivative ∂x/∂φ (latitude).
pub dx_dphi: f64,
/// Partial derivative ∂y/∂λ (longitude).
pub dy_dlam: f64,
/// Partial derivative ∂y/∂φ (latitude).
pub dy_dphi: f64,
/// `true` when computed via exact AD; `false` if a finite-difference
/// fallback was used instead.
pub is_exact: bool,
}
impl FactorsExact {
/// Compute distortion factors from the 2×2 Jacobian and geodetic latitude.
///
/// The four partial derivatives are obtained by evaluating a
/// `ProjectGeneric::project_fwd_generic` with `Dual1<2>` variables and
/// reading the `.d` arrays:
///
/// ```text
/// lam_d = Dual1::<2>::variable(lam, 0);
/// phi_d = Dual1::<2>::variable(phi, 1);
/// let (x, y) = proj.project_fwd_generic(lam_d, phi_d)?;
/// // dx_dlam = x.d[0], dx_dphi = x.d[1]
/// // dy_dlam = y.d[0], dy_dphi = y.d[1]
/// ```
///
/// `a` is the semi-major axis in the projection's working units (usually
/// normalised to 1.0 inside the engine, or `6_378_137.0` m for WGS 84).
pub fn from_jacobian(
dx_dlam: f64,
dx_dphi: f64,
dy_dlam: f64,
dy_dphi: f64,
phi: f64,
a: f64,
) -> Self {
let cos_phi = phi.cos().abs().max(1e-10);
// Meridional scale: h = ‖(dx_dphi, dy_dphi)‖ / a
let h = (dx_dphi * dx_dphi + dy_dphi * dy_dphi).sqrt() / a;
// Parallel scale: k = ‖(dx_dlam, dy_dlam)‖ / (a·cos φ)
let k = (dx_dlam * dx_dlam + dy_dlam * dy_dlam).sqrt() / (a * cos_phi);
// Areal scale: s = |det J| / (a²·cos φ)
let det = dx_dlam * dy_dphi - dx_dphi * dy_dlam;
let s = det.abs() / (a * a * cos_phi);
// Angular distortion: ω = 2·arcsin(|h−k|/(h+k))
let omega = if (h + k) > 1e-14 {
2.0 * ((h - k).abs() / (h + k)).min(1.0).asin()
} else {
0.0
};
// Meridian convergence: angle between grid north and true north
let conv = dx_dphi.atan2(dy_dphi);
Self {
meridian_scale: h,
parallel_scale: k,
areal_scale: s,
angular_distortion: omega,
meridian_convergence: conv,
dx_dlam,
dx_dphi,
dy_dlam,
dy_dphi,
is_exact: true,
}
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn spherical_mercator_at_equator() {
// At phi=0: h=k=1, s=1, omega=0, conv=0
// Spherical Mercator Jacobian at (lam, phi=0):
// x = lam => dx/dlam=1, dx/dphi=0
// y = ln(tan(pi/4 + phi/2)) => dy/dlam=0, dy/dphi=sec(0)=1
let fe = FactorsExact::from_jacobian(1.0, 0.0, 0.0, 1.0, 0.0, 1.0);
assert!(
(fe.meridian_scale - 1.0).abs() < 1e-12,
"h={}",
fe.meridian_scale
);
assert!(
(fe.parallel_scale - 1.0).abs() < 1e-12,
"k={}",
fe.parallel_scale
);
assert!((fe.areal_scale - 1.0).abs() < 1e-12, "s={}", fe.areal_scale);
assert!(
fe.angular_distortion.abs() < 1e-12,
"omega={}",
fe.angular_distortion
);
assert!(fe.is_exact);
}
#[test]
fn spherical_mercator_at_phi30() {
use core::f64::consts::FRAC_PI_6;
let phi = FRAC_PI_6; // 30 deg
let sec_phi = 1.0 / phi.cos();
// Spherical Mercator Jacobian at phi=30deg (a=1):
// dx/dlam=1, dx/dphi=0, dy/dlam=0, dy/dphi=sec(phi)
// h = sec(phi), k = 1/cos(phi) = sec(phi) (conformal => h=k)
// s = |det J| / (a^2 * cos phi) = sec(phi) / cos(phi) = sec^2(phi)
// Mercator is conformal but NOT equal-area; s != 1.
let fe = FactorsExact::from_jacobian(1.0, 0.0, 0.0, sec_phi, phi, 1.0);
assert!(
(fe.meridian_scale - sec_phi).abs() < 1e-12,
"h={}",
fe.meridian_scale
);
assert!(
(fe.parallel_scale - sec_phi).abs() < 1e-12,
"k={}",
fe.parallel_scale
);
let expected_s = sec_phi * sec_phi;
assert!(
(fe.areal_scale - expected_s).abs() < 1e-12,
"s={}",
fe.areal_scale
);
assert!(
fe.angular_distortion.abs() < 1e-12,
"omega={}",
fe.angular_distortion
);
}
}