geodesy/inner_op/
gridshift.rs

1/// Datum shift using grid interpolation.
2use crate::authoring::*;
3
4// ----- F O R W A R D --------------------------------------------------------------
5
6fn fwd(op: &Op, ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
7    let grids = &op.params.grids;
8    let use_null_grid = op.params.boolean("null_grid");
9
10    let mut successes = 0_usize;
11    let n = operands.len();
12
13    // Nothing to do?
14    if grids.is_empty() {
15        return n;
16    }
17
18    let ctx = Some(ctx);
19
20    for i in 0..n {
21        let mut coord = operands.get_coord(i);
22
23        if let Some(d) = grids_at(ctx, grids, coord, use_null_grid) {
24            // Geoid
25            if grids[0].bands() == 1 {
26                coord[2] -= d[0];
27                operands.set_coord(i, &coord);
28                successes += 1;
29
30                continue;
31            }
32
33            // Datum shift
34            let d = d.arcsec_to_radians();
35            coord[0] += d[0];
36            coord[1] += d[1];
37            operands.set_coord(i, &coord);
38            successes += 1;
39
40            continue;
41        }
42
43        // No grid contained the point, so we stomp on the coordinate
44        operands.set_coord(i, &Coor4D::nan());
45    }
46
47    successes
48}
49
50// ----- I N V E R S E --------------------------------------------------------------
51
52fn inv(op: &Op, ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
53    let grids = &op.params.grids;
54    let use_null_grid = op.params.boolean("null_grid");
55
56    let mut successes = 0_usize;
57    let n = operands.len();
58
59    // Nothing to do?
60    if grids.is_empty() {
61        return n;
62    }
63
64    let ctx = Some(ctx);
65
66    'points: for i in 0..n {
67        let mut coord = operands.get_coord(i);
68        if let Some(t) = grids_at(ctx, grids, coord, use_null_grid) {
69            // Geoid
70            if grids[0].bands() == 1 {
71                coord[2] += t[0];
72                operands.set_coord(i, &coord);
73                successes += 1;
74                continue;
75            }
76
77            // Inverse case datum shift - iteration needed
78            let t = t.arcsec_to_radians();
79            let mut t = coord - t;
80            for _ in 0..10 {
81                if let Some(t2) = grids_at(ctx, grids, t, use_null_grid) {
82                    let d = t - coord + t2.arcsec_to_radians();
83                    t = t - d;
84                    if d[0].hypot(d[1]) < 1e-12 {
85                        operands.set_coord(i, &t);
86                        successes += 1;
87                        continue 'points;
88                    }
89                    continue;
90                }
91
92                // The iteration has wandered off the grids, so we stomp
93                // on the coordinate and go on with the next
94                operands.set_coord(i, &Coor4D::nan());
95                continue 'points;
96            }
97        }
98    }
99
100    successes
101}
102
103// ----- C O N S T R U C T O R ------------------------------------------------------
104
105#[rustfmt::skip]
106pub const GAMUT: [OpParameter; 3] = [
107    OpParameter::Flag { key: "inv" },
108    OpParameter::Texts { key: "grids", default: None },
109    OpParameter::Real { key: "padding", default: Some(0.5) },
110];
111
112pub fn new(parameters: &RawParameters, ctx: &dyn Context) -> Result<Op, Error> {
113    let def = &parameters.instantiated_as;
114    let mut params = ParsedParameters::new(parameters, &GAMUT)?;
115
116    for mut grid_name in params.texts("grids")?.clone() {
117        let optional = grid_name.starts_with('@');
118        if optional {
119            grid_name = grid_name.trim_start_matches('@').to_string();
120        }
121
122        if grid_name == "null" {
123            params.boolean.insert("null_grid");
124            break; // ignore any additional grids after a null grid
125        }
126
127        match ctx.get_grid(&grid_name) {
128            Ok(grid) => params.grids.push(grid),
129            Err(e) => {
130                if !optional {
131                    return Err(e);
132                }
133            }
134        }
135    }
136
137    let fwd = InnerOp(fwd);
138    let inv = InnerOp(inv);
139    let descriptor = OpDescriptor::new(def, fwd, Some(inv));
140
141    Ok(Op {
142        descriptor,
143        params,
144        steps: None,
145    })
146}
147
148// ----- T E S T S ------------------------------------------------------------------
149
150//#[cfg(with_plain)]
151#[cfg(test)]
152mod tests {
153    use super::*;
154    use crate::coordinate::AngularUnits;
155
156    #[test]
157    fn gridshift() -> Result<(), Error> {
158        let mut ctx = Plain::default();
159        let op = ctx.op("gridshift grids=test.datum")?;
160        let cph = Coor4D::geo(55., 12., 0., 0.);
161        let mut data = [cph];
162
163        ctx.apply(op, Fwd, &mut data)?;
164        let res = data[0].to_geo();
165        assert!((res[0] - 55.015278).abs() < 1e-6);
166        assert!((res[1] - 12.003333).abs() < 1e-6);
167
168        ctx.apply(op, Inv, &mut data)?;
169        assert!((data[0][0] - cph[0]).abs() < 1e-10);
170        assert!((data[0][1] - cph[1]).abs() < 1e-10);
171
172        Ok(())
173    }
174
175    #[test]
176    fn ntv2() -> Result<(), Error> {
177        let mut ctx = Plain::default();
178        let op = ctx.op("gridshift grids=100800401.gsb")?;
179        let bcn = Coor2D::geo(41.3874, 2.1686);
180        let mut data = [bcn];
181
182        ctx.apply(op, Fwd, &mut data)?;
183        let res = data[0].to_geo();
184        assert!((res[0] - 41.38627500250805).abs() < 1e-8);
185        assert!((res[1] - 2.167450821894838).abs() < 1e-8);
186
187        ctx.apply(op, Inv, &mut data)?;
188        assert!((data[0][0] - bcn[0]).abs() < 1e-10);
189        assert!((data[0][1] - bcn[1]).abs() < 1e-10);
190
191        Ok(())
192    }
193
194    #[test]
195    fn multiple_grids() -> Result<(), Error> {
196        let mut ctx = Plain::default();
197        let op = ctx.op("gridshift grids=test.datum, test.datum")?;
198        let cph = Coor4D::geo(55., 12., 0., 0.);
199        let mut data = [cph];
200
201        ctx.apply(op, Fwd, &mut data)?;
202        let res = data[0].to_geo();
203        assert!((res[0] - 55.015278).abs() < 1e-6);
204        assert!((res[1] - 12.003333).abs() < 1e-6);
205
206        // Check that the reference counting works as expected
207        Plain::clear_grids();
208        Plain::clear_grids();
209
210        ctx.apply(op, Inv, &mut data)?;
211        assert!((data[0][0] - cph[0]).abs() < 1e-10);
212        assert!((data[0][1] - cph[1]).abs() < 1e-10);
213
214        Ok(())
215    }
216
217    #[test]
218    fn fails_without_null_grid() -> Result<(), Error> {
219        let mut ctx = Plain::default();
220        let op = ctx.op("gridshift grids=test.datum")?;
221
222        let ldn = Coor4D::geo(51.505, -0.09, 0., 0.);
223        let mut data = [ldn];
224
225        let successes = ctx.apply(op, Fwd, &mut data)?;
226        assert_eq!(successes, 0);
227        assert!(data[0][0].is_nan());
228        assert!(data[0][1].is_nan());
229
230        Ok(())
231    }
232
233    #[test]
234    fn passes_with_null_grid() -> Result<(), Error> {
235        let mut ctx = Plain::default();
236        let op = ctx.op("gridshift grids=test.datum, @null")?;
237
238        let ldn = Coor4D::geo(51.505, -0.09, 0., 0.);
239        let mut data = [ldn];
240
241        let successes = ctx.apply(op, Fwd, &mut data)?;
242        let res = data[0].to_geo();
243        assert_eq!(successes, 1);
244        assert_eq!(res[0], 51.505);
245        assert_eq!(res[1], -0.09);
246
247        let successes = ctx.apply(op, Inv, &mut data)?;
248        assert_eq!(successes, 1);
249        assert!((data[0][0] - ldn[0]).abs() < 1e-10);
250        assert!((data[0][1] - ldn[1]).abs() < 1e-10);
251
252        Ok(())
253    }
254
255    #[test]
256    fn optional_grid() -> Result<(), Error> {
257        let mut ctx = Plain::default();
258        let op = ctx.op("gridshift grids=@test_subset.datum, @missing.gsb, test.datum")?;
259
260        // Check that the Arc reference counting and the clear-GRIDS-by-instantiation
261        // works together correctly (i.e. the Arcs used by op are kept alive when
262        // GRIDS is cleared on the instantiation of ctx2)
263        let mut ctx2 = Plain::default();
264        let op2 = ctx2.op("gridshift grids=@test_subset.datum, @missing.gsb, test.datum")?;
265
266        // Copenhagen is outside of the (optional, but present, subset grid)
267        let cph = Coor2D::geo(55., 12.);
268        let mut data = [cph];
269        let mut data2 = [cph];
270
271        ctx.apply(op, Fwd, &mut data)?;
272        let res = data[0].to_geo();
273        assert!((res[0] - 55.015278).abs() < 1e-6);
274        assert!((res[1] - 12.003333).abs() < 1e-6);
275
276        // Same procedure on the other context gives same result
277        ctx2.apply(op2, Fwd, &mut data2)?;
278        assert_eq!(data, data2);
279
280        ctx.apply(op, Inv, &mut data)?;
281        assert!((data[0][0] - cph[0]).abs() < 1e-10);
282        assert!((data[0][1] - cph[1]).abs() < 1e-10);
283
284        // Same procedure on the other context gives same result
285        ctx2.apply(op2, Inv, &mut data2)?;
286        assert_eq!(data, data2);
287
288        // Havnebyen (a small town with a large geodetic installation) is inside the subset grid
289        let haby = Coor4D::geo(55.97, 11.33, 0., 0.);
290        let mut data = [haby];
291        let expected_correction = Coor4D([11.331, 55.971, 0., 0.]);
292        ctx.apply(op, Fwd, &mut data)?;
293        let correction = ((data[0] - haby) * Coor4D([3600., 3600., 3600., 3600.])).to_degrees();
294        assert!((correction - expected_correction)[0].abs() < 1e-6);
295        assert!((correction - expected_correction)[1].abs() < 1e-6);
296
297        Ok(())
298    }
299
300    #[test]
301    fn missing_grid() -> Result<(), Error> {
302        let mut ctx = Plain::default();
303        let op = ctx.op("gridshift grids=missing.gsb");
304        assert!(op.is_err());
305
306        Ok(())
307    }
308}
309
310// See additional tests in src/grid/mod.rs
geodesy/inner_op/gridshift.rs

geodesy/inner_op/
gridshift.rs
