Skip to content

Commit 7f1ec6a

Browse files
busstoptaktikbusstoptaktik
authored
Merge pull request #109 from busstoptaktik/work
- Streamlined CoorND operator macros - All projections now use xy()/set_xy() methods
2 parents b2a1560 + 0f3fcd0 commit 7f1ec6a

File tree

10 files changed

+101
-166
lines changed

10 files changed

+101
-166
lines changed

‎src/coordinate/coor2d.rs

+1-1
Original file line numberDiff line numberDiff line change
@@ -60,7 +60,7 @@ impl Coor2D {
6060
pub fn iso_dm(latitude: f64, longitude: f64) -> Coor2D {
6161
let longitude = angular::iso_dm_to_dd(longitude);
6262
let latitude = angular::iso_dm_to_dd(latitude);
63-
Coor2D([longitude.to_radians(), latitude.to_radians()])
63+
Coor2D::geo(latitude, longitude)
6464
}
6565

6666
/// A `Coor2D` from latitude/longitude in the ISO-6709 DDDMMSS.sssss format.

‎src/coordinate/coor32.rs

+1-1
Original file line numberDiff line numberDiff line change
@@ -66,7 +66,7 @@ impl Coor32 {
6666
pub fn iso_dm(latitude: f64, longitude: f64) -> Coor32 {
6767
let longitude = angular::iso_dm_to_dd(longitude);
6868
let latitude = angular::iso_dm_to_dd(latitude);
69-
Coor32([longitude.to_radians() as f32, latitude.to_radians() as f32])
69+
Coor32::geo(latitude, longitude)
7070
}
7171

7272
/// A `Coor32` from latitude/longitude/height/time, with the

‎src/coordinate/coor3d.rs

+1-1
Original file line numberDiff line numberDiff line change
@@ -61,7 +61,7 @@ impl Coor3D {
6161
pub fn iso_dm(latitude: f64, longitude: f64, height: f64) -> Coor3D {
6262
let longitude = angular::iso_dm_to_dd(longitude);
6363
let latitude = angular::iso_dm_to_dd(latitude);
64-
Coor3D([longitude.to_radians(), latitude.to_radians(), height])
64+
Coor3D::geo(latitude, longitude, height)
6565
}
6666

6767
/// A `Coor3D` from latitude/longitude/height, with

‎src/coordinate/coor4d.rs

+1-1
Original file line numberDiff line numberDiff line change
@@ -62,7 +62,7 @@ impl Coor4D {
6262
pub fn iso_dm(latitude: f64, longitude: f64, height: f64, time: f64) -> Coor4D {
6363
let longitude = angular::iso_dm_to_dd(longitude);
6464
let latitude = angular::iso_dm_to_dd(latitude);
65-
Coor4D([longitude.to_radians(), latitude.to_radians(), height, time])
65+
Coor4D::geo(latitude, longitude, height, time)
6666
}
6767

6868
/// A `Coor4D` from latitude/longitude/height/time, with the

‎src/coordinate/tuple.rs

+21-49
Original file line numberDiff line numberDiff line change
@@ -85,55 +85,27 @@ macro_rules! coord_operator {
8585
};
8686
}
8787

88-
coord_operator!(Coor4D, Coor4D, coor4d, Add, +, add);
89-
coord_operator!(Coor4D, Coor4D, coor4d, Sub, -, sub);
90-
coord_operator!(Coor4D, Coor4D, coor4d, Mul, *, mul);
91-
coord_operator!(Coor4D, Coor4D, coor4d, Div, /, div);
92-
93-
coord_operator!(Coor4D, &Coor4D, coor4d, Add, +, add);
94-
coord_operator!(Coor4D, &Coor4D, coor4d, Sub, -, sub);
95-
coord_operator!(Coor4D, &Coor4D, coor4d, Mul, *, mul);
96-
coord_operator!(Coor4D, &Coor4D, coor4d, Div, /, div);
97-
98-
coord_operator!(Coor3D, Coor3D, coor3d, Add, +, add);
99-
coord_operator!(Coor3D, Coor3D, coor3d, Sub, -, sub);
100-
coord_operator!(Coor3D, Coor3D, coor3d, Mul, *, mul);
101-
coord_operator!(Coor3D, Coor3D, coor3d, Div, /, div);
102-
103-
coord_operator!(Coor3D, &Coor3D, coor3d, Add, +, add);
104-
coord_operator!(Coor3D, &Coor3D, coor3d, Sub, -, sub);
105-
coord_operator!(Coor3D, &Coor3D, coor3d, Mul, *, mul);
106-
coord_operator!(Coor3D, &Coor3D, coor3d, Div, /, div);
107-
108-
coord_operator!(Coor2D, Coor2D, coor2d, Add, +, add);
109-
coord_operator!(Coor2D, Coor2D, coor2d, Sub, -, sub);
110-
coord_operator!(Coor2D, Coor2D, coor2d, Mul, *, mul);
111-
coord_operator!(Coor2D, Coor2D, coor2d, Div, /, div);
112-
113-
coord_operator!(Coor2D, &Coor2D, coor2d, Add, +, add);
114-
coord_operator!(Coor2D, &Coor2D, coor2d, Sub, -, sub);
115-
coord_operator!(Coor2D, &Coor2D, coor2d, Mul, *, mul);
116-
coord_operator!(Coor2D, &Coor2D, coor2d, Div, /, div);
117-
118-
coord_operator!(Coor2D, Coor32, coor2d, Add, +, add);
119-
coord_operator!(Coor2D, Coor32, coor2d, Sub, -, sub);
120-
coord_operator!(Coor2D, Coor32, coor2d, Mul, *, mul);
121-
coord_operator!(Coor2D, Coor32, coor2d, Div, /, div);
122-
123-
coord_operator!(Coor2D, &Coor32, coor2d, Add, +, add);
124-
coord_operator!(Coor2D, &Coor32, coor2d, Sub, -, sub);
125-
coord_operator!(Coor2D, &Coor32, coor2d, Mul, *, mul);
126-
coord_operator!(Coor2D, &Coor32, coor2d, Div, /, div);
127-
128-
coord_operator!(Coor32, Coor32, coor32, Add, +, add);
129-
coord_operator!(Coor32, Coor32, coor32, Sub, -, sub);
130-
coord_operator!(Coor32, Coor32, coor32, Mul, *, mul);
131-
coord_operator!(Coor32, Coor32, coor32, Div, /, div);
132-
133-
coord_operator!(Coor32, &Coor32, coor32, Add, +, add);
134-
coord_operator!(Coor32, &Coor32, coor32, Sub, -, sub);
135-
coord_operator!(Coor32, &Coor32, coor32, Mul, *, mul);
136-
coord_operator!(Coor32, &Coor32, coor32, Div, /, div);
88+
// Generate the vector space operators Add, Sub, Mul, Div for $type
89+
macro_rules! all_coord_operators {
90+
($type:ty, $othertype:ty, $typemacro:ident) => {
91+
coord_operator!($type, $othertype, $typemacro, Add, +, add);
92+
coord_operator!($type, $othertype, $typemacro, Sub, -, sub);
93+
coord_operator!($type, $othertype, $typemacro, Mul, *, mul);
94+
coord_operator!($type, $othertype, $typemacro, Div, /, div);
95+
};
96+
}
97+
98+
all_coord_operators!(Coor4D, &Coor4D, coor4d);
99+
all_coord_operators!(Coor3D, &Coor3D, coor3d);
100+
all_coord_operators!(Coor2D, &Coor2D, coor2d);
101+
all_coord_operators!(Coor2D, &Coor32, coor2d);
102+
all_coord_operators!(Coor32, &Coor32, coor32);
103+
104+
all_coord_operators!(Coor4D, Coor4D, coor4d);
105+
all_coord_operators!(Coor3D, Coor3D, coor3d);
106+
all_coord_operators!(Coor2D, Coor2D, coor2d);
107+
all_coord_operators!(Coor2D, Coor32, coor2d);
108+
all_coord_operators!(Coor32, Coor32, coor32);
137109

138110
/// CoordinateTuple is the ISO-19111 atomic spatial/spatiotemporal
139111
/// referencing element. So loosely speaking, a CoordinateSet is a

‎src/inner_op/laea.rs

+31-48
Original file line numberDiff line numberDiff line change
@@ -5,8 +5,6 @@ use crate::authoring::*;
55
use std::f64::consts::FRAC_PI_2;
66
const EPS10: f64 = 1e-10;
77

8-
// ----- C O M M O N -------------------------------------------------------------------
9-
108
// ----- F O R W A R D -----------------------------------------------------------------
119

1210
fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
@@ -40,49 +38,41 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
4038

4139
// The polar aspects are fairly simple
4240
if north_polar || south_polar {
41+
let sign = if north_polar { -1.0 } else { 1.0 };
4342
for i in 0..n {
44-
let mut coord = operands.get_coord(i);
45-
let sign = if north_polar { -1.0 } else { 1.0 };
46-
47-
let lat = coord[1];
48-
let lon = coord[0];
43+
let (lon, lat) = operands.xy(i);
4944
let (sin_lon, cos_lon) = (lon - lon_0).sin_cos();
5045

5146
let q = ancillary::qs(lat.sin(), e);
5247
let rho = a * (qp + sign * q).sqrt();
5348

54-
coord[0] = x_0 + rho * sin_lon;
55-
coord[1] = y_0 + sign * rho * cos_lon;
56-
operands.set_coord(i, &coord);
49+
let easting = x_0 + rho * sin_lon;
50+
let northing = y_0 + sign * rho * cos_lon;
51+
operands.set_xy(i, easting, northing);
5752
successes += 1;
5853
}
5954
return successes;
6055
}
6156

57+
// Either equatorial or oblique aspects
6258
for i in 0..n {
63-
let mut coord = operands.get_coord(i);
64-
let lon = coord[0];
65-
let lat = coord[1];
59+
let (lon, lat) = operands.xy(i);
6660
let (sin_lon, cos_lon) = (lon - lon_0).sin_cos();
6761

6862
// Authalic latitude, 𝜉
6963
let xi = (ancillary::qs(lat.sin(), e) / qp).asin();
70-
7164
let (sin_xi, cos_xi) = xi.sin_cos();
65+
7266
let b = if oblique {
7367
let factor = 1.0 + sin_xi_0 * sin_xi + (cos_xi_0 * cos_xi * cos_lon);
7468
rq * (2.0 / factor).sqrt()
7569
} else {
7670
1.0
7771
};
7872

79-
// Easting
80-
coord[0] = x_0 + (b * d) * (cos_xi * sin_lon);
81-
82-
// Northing
83-
coord[1] = y_0 + (b / d) * (cos_xi_0 * sin_xi - sin_xi_0 * cos_xi * cos_lon);
84-
operands.set_coord(i, &coord);
85-
73+
let easting = x_0 + (b * d) * (cos_xi * sin_lon);
74+
let northing = y_0 + (b / d) * (cos_xi_0 * sin_xi - sin_xi_0 * cos_xi * cos_lon);
75+
operands.set_xy(i, easting, northing);
8676
successes += 1;
8777
}
8878

@@ -122,39 +112,33 @@ fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
122112
let mut successes = 0_usize;
123113
let n = operands.len();
124114

125-
// The polar aspects are not as simple as in the forward case
115+
// The polar aspects are not quite as simple as in the forward case
126116
if north_polar || south_polar {
117+
let sign = if north_polar { -1.0 } else { 1.0 };
127118
for i in 0..n {
128-
let mut coord = operands.get_coord(i);
129-
let sign = if north_polar { -1.0 } else { 1.0 };
130-
131-
let x = coord[0];
132-
let y = coord[1];
119+
let (x, y) = operands.xy(i);
133120
let rho = (x - x_0).hypot(y - y_0);
134121

135122
// The authalic latitude is a bit convoluted
136123
let denom = a * a * (1.0 - ((1.0 - es) / (2.0 * e)) * ((1.0 - e) / (1.0 + e)).ln());
137124
let xi = (-sign) * (1.0 - rho * rho / denom);
138125

139-
coord[0] = lon_0 + (x - x_0).atan2(sign * (y - y_0));
140-
coord[1] = ellps.latitude_authalic_to_geographic(xi, &authalic);
141-
operands.set_coord(i, &coord);
126+
let lon = lon_0 + (x - x_0).atan2(sign * (y - y_0));
127+
let lat = ellps.latitude_authalic_to_geographic(xi, &authalic);
128+
operands.set_xy(i, lon, lat);
142129
successes += 1;
143130
}
144131
return successes;
145132
}
146133

134+
// Either equatorial or oblique aspects
147135
for i in 0..n {
148-
let mut coord = operands.get_coord(i);
149-
let x = coord[0];
150-
let y = coord[1];
136+
let (x, y) = operands.xy(i);
151137
let rho = ((x - x_0) / d).hypot(d * (y - y_0));
152138

153139
// A bit of reality hardening ported from the PROJ implementation
154140
if rho < EPS10 {
155-
coord[0] = lon_0;
156-
coord[1] = lat_0;
157-
operands.set_coord(i, &coord);
141+
operands.set_xy(i, lon_0, lat_0);
158142
successes += 1;
159143
continue;
160144
}
@@ -163,23 +147,22 @@ fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
163147
let asin_argument = 0.5 * rho / rq;
164148
if asin_argument.abs() > 1.0 {
165149
debug!("LAEA: ({x}, {y}) outside domain");
166-
coord[0] = f64::NAN;
167-
coord[1] = f64::NAN;
168-
operands.set_coord(i, &coord);
150+
operands.set_xy(i, f64::NAN, f64::NAN);
169151
continue;
170152
}
171153

172154
let c = 2.0 * asin_argument.asin();
173155
let (sin_c, cos_c) = c.sin_cos();
156+
174157
// The authalic latitude, 𝜉
175158
let xi = (cos_c * sin_xi_0 + (d * (y - y_0) * sin_c * cos_xi_0) / rho).asin();
176-
coord[1] = ellps.latitude_authalic_to_geographic(xi, &authalic);
159+
160+
let lat = ellps.latitude_authalic_to_geographic(xi, &authalic);
177161

178162
let num = (x - x_0) * sin_c;
179163
let denom = d * rho * cos_xi_0 * cos_c - d * d * (y - y_0) * sin_xi_0 * sin_c;
180-
coord[0] = num.atan2(denom) + lon_0;
181-
operands.set_coord(i, &coord);
182-
164+
let lon = num.atan2(denom) + lon_0;
165+
operands.set_xy(i, lon, lat);
183166
successes += 1;
184167
}
185168

@@ -219,12 +202,12 @@ pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result<Op, Error>
219202

220203
let polar = (t - FRAC_PI_2).abs() < EPS10;
221204
let north = polar && (t > 0.0);
222-
let equatoreal = !polar && t < EPS10;
223-
let oblique = !polar && !equatoreal;
224-
match (polar, equatoreal, north) {
205+
let equatorial = !polar && t < EPS10;
206+
let oblique = !polar && !equatorial;
207+
match (polar, equatorial, north) {
225208
(true, _, true) => params.boolean.insert("north_polar"),
226209
(true, _, false) => params.boolean.insert("south_polar"),
227-
(_, true, _) => params.boolean.insert("equatoreal"),
210+
(_, true, _) => params.boolean.insert("equatorial"),
228211
_ => params.boolean.insert("oblique"),
229212
};
230213

@@ -247,7 +230,7 @@ pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result<Op, Error>
247230
// D in the IOGP text
248231
let d = if oblique {
249232
a * (cos_phi_0 / (1.0 - es * sin_phi_0 * sin_phi_0).sqrt()) / (rq * xi_0.cos())
250-
} else if equatoreal {
233+
} else if equatorial {
251234
rq.recip()
252235
} else {
253236
a

‎src/inner_op/lcc.rs

+19-24
Original file line numberDiff line numberDiff line change
@@ -25,25 +25,23 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
2525
let length = operands.len();
2626

2727
for i in 0..length {
28-
let mut coord = operands.get_coord(i);
29-
let lam = coord[0] - lon_0;
30-
let phi = coord[1];
28+
let (mut lam, phi) = operands.xy(i);
29+
lam -= lon_0;
3130
let mut rho = 0.;
3231

3332
// Close to one of the poles?
3433
if (phi.abs() - FRAC_PI_2).abs() < EPS10 {
3534
if phi * n <= 0. {
36-
coord = Coor4D::nan();
37-
operands.set_coord(i, &coord);
35+
operands.set_coord(i, &Coor4D::nan());
3836
continue;
3937
}
4038
} else {
4139
rho = c * crate::math::ancillary::ts(phi.sin_cos(), e).powf(n);
4240
}
4341
let sc = (lam * n).sin_cos();
44-
coord[0] = a * k_0 * rho * sc.0 + x_0;
45-
coord[1] = a * k_0 * (rho0 - rho * sc.1) + y_0;
46-
operands.set_coord(i, &coord);
42+
let x = a * k_0 * rho * sc.0 + x_0;
43+
let y = a * k_0 * (rho0 - rho * sc.1) + y_0;
44+
operands.set_xy(i, x, y);
4745
successes += 1;
4846
}
4947
successes
@@ -64,20 +62,19 @@ fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
6462
return 0;
6563
};
6664
let mut successes = 0_usize;
67-
let length = operands.len();
6865

69-
for i in 0..length {
70-
let mut coord = operands.get_coord(i);
71-
let mut x = (coord[0] - x_0) / (a * k_0);
72-
let mut y = rho0 - (coord[1] - y_0) / (a * k_0);
66+
for i in 0..operands.len() {
67+
let (mut x, mut y) = operands.xy(i);
68+
x = (x - x_0) / (a * k_0);
69+
y = rho0 - (y - y_0) / (a * k_0);
7370

7471
let mut rho = x.hypot(y);
7572

76-
// On one of the poles
73+
// On one of the poles?
7774
if rho == 0. {
78-
coord[0] = 0.;
79-
coord[1] = FRAC_PI_2.copysign(n);
80-
operands.set_coord(i, &coord);
75+
let lon = 0.;
76+
let lat = FRAC_PI_2.copysign(n);
77+
operands.set_xy(i, lon, lat);
8178
successes += 1;
8279
continue;
8380
}
@@ -90,15 +87,13 @@ fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
9087
}
9188

9289
let ts0 = (rho / c).powf(1. / n);
93-
let phi = crate::math::ancillary::pj_phi2(ts0, e);
94-
if phi.is_infinite() || phi.is_nan() {
95-
coord = Coor4D::nan();
96-
operands.set_coord(i, &coord);
90+
let lat = crate::math::ancillary::pj_phi2(ts0, e);
91+
if lat.is_infinite() || lat.is_nan() {
92+
operands.set_coord(i, &Coor4D::nan());
9793
continue;
9894
}
99-
coord[0] = x.atan2(y) / n + lon_0;
100-
coord[1] = phi;
101-
operands.set_coord(i, &coord);
95+
let lon = x.atan2(y) / n + lon_0;
96+
operands.set_xy(i, lon, lat);
10297
successes += 1;
10398
}
10499
successes

‎src/inner_op/merc.rs

+2-4
Original file line numberDiff line numberDiff line change
@@ -13,8 +13,7 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
1313
let lon_0 = op.params.lon(0);
1414

1515
let mut successes = 0_usize;
16-
let length = operands.len();
17-
for i in 0..length {
16+
for i in 0..operands.len() {
1817
let (lon, lat) = operands.xy(i);
1918

2019
let easting = (lon - lon_0) * k_0 * a - x_0;
@@ -40,8 +39,7 @@ fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
4039
let lon_0 = op.params.lon(0);
4140

4241
let mut successes = 0_usize;
43-
let length = operands.len();
44-
for i in 0..length {
42+
for i in 0..operands.len() {
4543
let (mut x, mut y) = operands.xy(i);
4644

4745
// Easting -> Longitude

0 commit comments

Comments
 (0)