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
159
160
use crate::prelude::*;
/// The Jacobian matrix for investigation of the geometrical properties of
/// map projections. Can be converted to the more easily digestible
/// [`Factors`] struct, using the `.factors()` method.
#[derive(Debug, Default)]
#[rustfmt::skip]
pub struct Jacobian {
// The geographical coordinates (in degrees) of the evaluation point
pub latitude: f64,
pub longitude: f64,
// The derivatives at the evaluation point of easting (x) and northing (y)
// with respect to longitude (lam) and latitude (phi)
pub dx_dlam: f64,
pub dy_dlam: f64,
pub dx_dphi: f64,
pub dy_dphi: f64,
// The ellipsoid on which the jacobian is evaluated
pub ellps: Ellipsoid
}
/// Geometrical properties of map projections. The `Factors` struct is
/// derived from the more fundamental `Jacobian`, using its `.factors()`
/// method
#[derive(Debug, Default)]
#[rustfmt::skip]
pub struct Factors {
// Scalar factors // Common textbook designation
pub meridional_scale: f64, // h
pub parallel_scale: f64, // k
pub areal_scale: f64, // s
// Angular factors
pub angular_distortion: f64, // ω
pub meridian_parallel_angle: f64, // θ'
pub meridian_convergence: f64, // α
// Tissot indicatrix
pub tissot_semimajor: f64, // a
pub tissot_semiminor: f64, // b
}
impl Jacobian {
/// Compute the Jacobian matrix for the map projection represented by `op`.
/// The `scale` parameters define the scaling from input units to degrees,
/// and from output units to metres. Hence, if input is in radians and
/// output in feet, `scale`should be set to `[1f64.to_degrees(), 0.3048]`.
/// The `swap` parameters indicate whether input and output is swapped
/// with respect to the Rust Geodesy internal GIS convention of "longitude
/// before latitude, and easting before northing". Hence, if input is in
/// the geographical convention of latitude/longitude, and output is in
/// easting/northing, `swap` should be set to `[true, false]`. The `ellps`
/// parameter should be set to the relevant ellipsoid for the projection.
/// While it could be derived directly from the other input data in the
/// case of a single operation, the potential case of a pipeline operation
/// makes it necessary to actively select which one to use. In all but the
/// most demanding situations it is, however, probably fine to just select
/// `Ellipsoid::default()`, i.e. GRS80.
///
/// Mostly based on the PROJ function [pj_deriv](https://github.com/OSGeo/PROJ/blob/master/src/deriv.cpp),
#[rustfmt::skip]
pub fn new(ctx: &impl Context, op: OpHandle, scale: [f64; 2], swap: [bool; 2], ellps: Ellipsoid, at: Coor2D) -> Result<Jacobian, Error> {
// If we have input in degrees, we must multiply the output by a factor of 180/pi
// For user convenience, scale[0] is a "to degrees"-factor, i.e. scale[0]==1
// indicates degrees, whereas scale[0]==180/pi indicates that input angles are
// in radians.
// To convert input coordinates to radians, divide by `angular_scale`
let angular_scale = 1f64.to_degrees() / scale[0];
// If we have output in feet, we must multiply the output by a factor of 0.3048
// For user convenience, scale[1] is a "to metres"-factor, i.e. scale[1]==1
// indicates metres, whereas scale[1]==0.3048 indicates that output lengths
// are in feet, and scale[1]=201.168 indicates that output is in furlongs
let linear_scale = scale[1];
let h = 1e-5 * angular_scale;
let d = (4.0 * h * ellps.semimajor_axis()).recip() * linear_scale * angular_scale;
let mut coo = [Coor2D::origin(); 4];
let (e, n) = if swap[0] {(at[1], at[0])} else {(at[0], at[1])};
// Latitude and longitude in degrees for the return value
let latitude = n * scale[0];
let longitude = e * scale[0];
// North-east of POI
coo[0] = Coor2D::raw(e + h, n + h);
// South-east of POI
coo[1] = Coor2D::raw(e + h, n - h);
// South-west of POI
coo[2] = Coor2D::raw(e - h, n - h);
// North-west of POI
coo[3] = Coor2D::raw(e - h, n + h);
if swap[0] {
coo[0] = Coor2D::raw(coo[0][1], coo[0][0]);
coo[1] = Coor2D::raw(coo[1][1], coo[1][0]);
coo[2] = Coor2D::raw(coo[2][1], coo[2][0]);
coo[3] = Coor2D::raw(coo[3][1], coo[3][0]);
}
ctx.apply(op, Fwd, &mut coo)?;
// Handle output swapping
let (e, n) = if swap[1] {(1, 0)} else {(0, 1)};
// NE SE SW NW
let dx_dlam = (coo[0][e] + coo[1][e] - coo[2][e] - coo[3][e]) * d;
let dy_dlam = (coo[0][n] + coo[1][n] - coo[2][n] - coo[3][n]) * d;
let dx_dphi = (coo[0][e] - coo[1][e] - coo[2][e] + coo[3][e]) * d;
let dy_dphi = (coo[0][n] - coo[1][n] - coo[2][n] + coo[3][n]) * d;
Ok(Jacobian{latitude, longitude, dx_dlam, dy_dlam, dx_dphi, dy_dphi, ellps})
}
/// This closely follows the PROJ function pj_factors() and its friendly wrapper
/// proj_factors(), i.e. closely following Snyder's magnum opus
pub fn factors(&self) -> Factors {
let mut f = Factors::default();
let x_l = self.dx_dlam;
let y_l = self.dy_dlam;
let x_p = self.dx_dphi;
let y_p = self.dy_dphi;
let (s, c) = self.latitude.to_radians().sin_cos();
let es = self.ellps.eccentricity_squared();
// Linear scaling factors
let h = x_p.hypot(y_p);
let k = x_l.hypot(y_l) / c;
// Correction of linear scaling factors for ellipsoidal geometry
let t = 1. - es * s * s;
let n = t.sqrt();
let h = h * (t * n / (1. - es));
let k = k * n;
let r = t * t / (1. - es);
f.meridional_scale = h;
f.parallel_scale = k;
f.areal_scale = (y_p * x_l - x_p * y_l) * r / c;
// Tissot axes
let t = h * h + k * k;
let a = (t + 2. * f.areal_scale).sqrt();
let t = (t - 2. * f.areal_scale).clamp(0., f64::MAX).sqrt();
f.tissot_semiminor = 0.5 * (a - t);
f.tissot_semimajor = 0.5 * (a + t);
// Angular elements
f.meridian_parallel_angle = (f.areal_scale / (h * k)).clamp(-1., 1.).asin().to_degrees();
f.meridian_convergence = -x_p.atan2(y_p).to_degrees();
let a = f.tissot_semimajor;
let b = f.tissot_semiminor;
f.angular_distortion = 2. * ((a - b) / (a + b)).asin();
f
}
}