diff --git a/geodesy/datum/test_subset.datum b/geodesy/datum/test_subset.datum new file mode 100644 index 0000000..af19229 --- /dev/null +++ b/geodesy/datum/test_subset.datum @@ -0,0 +1,12 @@ +# Subset of test.datum, with slightly offset corrections +# for tests of multi-grid functionality +# +# e.g. at P: (56, 11) we should get the +# result (56.001, 11.001) + + +55.5 57.5 11. 13. 1. 1. + + 57.501 11.001 57.501 12.001 57.501 13.001 + 56.501 11.001 56.501 12.001 56.501 13.001 + 55.501 11.001 55.501 12.001 55.501 13.001 diff --git a/geodesy/deformation/another_test.deformation b/geodesy/deformation/another_test.deformation new file mode 100644 index 0000000..52f9694 --- /dev/null +++ b/geodesy/deformation/another_test.deformation @@ -0,0 +1,8 @@ +# Like test.deformation, but moved (10, 10) degrees to the north-east +64. 68. 18. 26. 1. 1. + + 68. 18. 0. 68. 19. 0. 68. 20. 0. 68. 21. 0. 68. 22. 0. 68. 23. 0. 68. 24. 0. 68. 25. 0. 68. 26. 0. + 67. 18. 0. 67. 19. 0. 67. 20. 0. 67. 21. 0. 67. 22. 0. 67. 23. 0. 67. 24. 0. 67. 25. 0. 67. 26. 0. + 66. 18. 0. 66. 19. 0. 66. 20. 0. 66. 21. 0. 66. 22. 0. 66. 23. 0. 66. 24. 0. 66. 25. 0. 66. 26. 0. + 65. 18. 0. 65. 19. 0. 65. 20. 0. 65. 21. 0. 65. 22. 0. 65. 23. 0. 65. 24. 0. 65. 25. 0. 65. 26. 0. + 64. 18. 0. 64. 19. 0. 64. 20. 0. 64. 21. 0. 64. 22. 0. 64. 23. 0. 64. 24. 0. 64. 25. 0. 64. 26. 0. diff --git a/ruminations/002-rumination.md b/ruminations/002-rumination.md index 3737e9b..3bbc228 100644 --- a/ruminations/002-rumination.md +++ b/ruminations/002-rumination.md @@ -3,9 +3,10 @@ ## Rumination 002: The missing manual Thomas Knudsen + Sean Rennie -2021-08-20. Last [revision](#document-history) 2023-10-19 +2021-08-20. Last [revision](#document-history) 2023-11-02 ### Abstract @@ -45,7 +46,7 @@ $ echo 553036. -124509 | kp "dms:in | geo:out" Architecturally, the operators in Rust Geodesy (`cart`, `tmerc`, `helmert` etc.) live below the API surface. This means they are not (and should not be) described in the API documentation over at [docs.rs](https://docs.rs/geodesy). Rather, their use should be documented in a separate *Rust Geodesy User's Guide*, a book which may materialize some day, as time permits, interest demands, and RG has matured and stabilized sufficiently. Until then, this *Rumination* will serve as stop gap for operator documentation. -A *Rust Geodesy Programmer's Guide* would probably also be useful, and wil definitely materialize before the next week with ten fridays. Until then, the [API documentation](https://docs.rs/geodesy), the [code examples](/examples), and the [architectural overview](/ruminations/000-rumination.md) may be useful. The RG transformation program `kp` is described in [RG Rumination 003](/ruminations/003-rumination.md). Its [source code](/src/bin/kp.rs) may also be of interest as study material for programmers. But since it is particularly useful for practical experimentation with RG operators, let's start with a *very* brief description of `kp`. +A *Rust Geodesy Programmer's Guide* would probably also be useful, and will definitely materialize before the next week with ten fridays. Until then, the [API documentation](https://docs.rs/geodesy), the [code examples](/examples), and the [architectural overview](/ruminations/000-rumination.md) may be useful. The RG transformation program `kp` is described in [RG Rumination 003](/ruminations/003-rumination.md). Its [source code](/src/bin/kp.rs) may also be of interest as study material for programmers. But since it is particularly useful for practical experimentation with RG operators, let's start with a *very* brief description of `kp`. ### A brief `kp` HOWTO @@ -304,7 +305,7 @@ The `gridshift` operator implements datum shifts by interpolation in correction | Parameter | Description | |-----------|-------------| | `inv` | Inverse operation: output-to-input datum. For 2-D and 3-D cases, this involves an iterative refinement, typically converging after less than 5 iterations | -| `grids` | Name of the grid file to use. RG supports only one file for each operation, but maintains the plural form of the `grids` option for alignment with the PROJ precedent | +| `grids` | Name of the grid files to use. RG supports multiple comma separated grids where the first one to contain the point is the one used. Grids are considered optional if they are prefixed with `@` and do not error the operator if they aren't available. Additionally the `@null` parameter can be specified as the last grid which will prevent errors in shifts from stomping on the coordinate. That is to say the coordinate passes through unchanged. | The `gridshift` operator has built in support for the **Gravsoft** grid format. Support for additional file formats depends on the `Context` in use. @@ -315,6 +316,10 @@ For grids with angular (geographical) spatial units, the corrections are suppose ```term geo:in | gridshift grids=ed50.datum | geo:out + +geo:in | gridshift grids=ed50.datum,@null | geo:out + +geo:in | gridshift grids=@not-available.gsb,ed50.datum | geo:out ``` **See also:** PROJ documentation, [`hgridshift`](https://proj.org/operations/transformations/hgridshift.html) and [`vgridshift`](https://proj.org/operations/transformations/vgridshift.html). RG combines the functionality of the two: The dimensionality of the grid determines whether a plane or a vertical transformation is carried out. @@ -650,7 +655,7 @@ Take a copy of one or more coordinate dimensions and push it onto the stack. If somerc lat_0=46.9524055555556 lon_0=7.43958333333333 k_0=1 x_0=2600000 y_0=1200000 ellps=bessel ``` -**See also:** [PROJ documentation](https://proj.org/operations/projections/somerc.html): _Swiss Oblique Mercator_. +**See also:** [PROJ documentation](https://proj.org/operations/projections/somerc.html): *Swiss Oblique Mercator*. Note: Rust Geodesy does not support modifying the ellipsoid with and `R` parameter, as PROJ does. @@ -738,4 +743,5 @@ Major revisions and additions: registered update on 2022-05-08. a large number of new operators have been included and described - 2023-07-09: dm and dms liberated from their NMEA overlord -- 2023-10-19: Add `somerc` operator description \ No newline at end of file +- 2023-10-19: Add `somerc` operator description +- 2023-11-02: Update `gridshift` operator description with multi, optional and null grid support diff --git a/src/inner_op/deformation.rs b/src/inner_op/deformation.rs index 54d4d60..8006c41 100644 --- a/src/inner_op/deformation.rs +++ b/src/inner_op/deformation.rs @@ -76,7 +76,7 @@ /// now introduce the designations *observed coordinates* for (X, Y, Z), and /// *canonical coordinates* for (X', Y', Z'). /// -/// What we wnat to do is to compute the canonical coordinates given the +/// What we want to do is to compute the canonical coordinates given the /// observed ones, by applying a correction based on the deformation grid. /// /// The deformation grid is georeferenced with respect to the *canonical system* @@ -143,12 +143,6 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { let deformation = rotate_and_integrate_velocity(v.scale(-1.), geo[0], geo[1], d); - // Outside of the grid? - stomp on the input coordinate and go on to the next - if v[0].is_nan() { - operands.set_coord(i, &Coor4D::nan()); - continue 'points; - } - // Finally apply the deformation to the input coordinate - or just // provide the raw correction if that was what was requested if raw { @@ -204,12 +198,6 @@ fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { let deformation = rotate_and_integrate_velocity(v, geo[0], geo[1], d); - // Outside of the grid? - stomp on the input coordinate and go on to the next - if v[0].is_nan() { - operands.set_coord(i, &Coor4D::nan()); - continue 'points; - } - // Finally apply the deformation to the input coordinate - or just // provide the raw correction if that was what was requested if raw { @@ -262,24 +250,34 @@ pub fn new(parameters: &RawParameters, ctx: &dyn Context) -> Result { )); } - for grid_name in params.texts("grids")?.clone() { - if grid_name.ends_with("@null") { + for mut grid_name in params.texts("grids")?.clone() { + let optional = grid_name.starts_with('@'); + if optional { + grid_name = grid_name.trim_start_matches('@').to_string(); + } + if grid_name == "null" { params.boolean.insert("null_grid"); - continue; + break; // ignore any additional grids after a null grid } + match ctx.get_grid(&grid_name) { + Ok(grid) => { + let n = grid.bands(); + if n != 3 { + return Err(Error::Unexpected { + message: "Bad dimensionality of deformation model grid".to_string(), + expected: "3".to_string(), + found: n.to_string(), + }); + } + params.grids.push(grid); + } - // TODO: Handle @optional grids - - let grid = ctx.get_grid(&grid_name)?; - let n = grid.bands(); - if n != 3 { - return Err(Error::Unexpected { - message: "Bad dimensionality of deformation model grid".to_string(), - expected: "3".to_string(), - found: n.to_string(), - }); + Err(e) => { + if !optional { + return Err(e); + } + } } - params.grids.push(grid); } let fwd = InnerOp(fwd); @@ -335,9 +333,12 @@ mod tests { let mut ctx = Plain::default(); let cph = Coor4D::geo(55., 12., 0., 0.); let test_deformation = include_str!("../../geodesy/deformation/test.deformation"); + let another_test_deformation = + include_str!("../../geodesy/deformation/another_test.deformation"); // Check that grid registration works ctx.register_resource("test.deformation", test_deformation); + ctx.register_resource("another_test.deformation", another_test_deformation); let buf = ctx.get_blob("test.deformation")?; let grid = BaseGrid::gravsoft(&buf)?; @@ -391,7 +392,8 @@ mod tests { assert!(cph.hypot3(&data[0]) < 1e-3); // Check the "raw" functionality - let op = ctx.op("deformation raw dt=1000 grids=test.deformation")?; + let op = + ctx.op("deformation raw dt=1000 grids=@another_test.deformation,test.deformation")?; // Forward direction let mut data = [cph]; @@ -408,6 +410,22 @@ mod tests { assert!((inv[3] - expected_length_of_correction) < 0.001); assert!((inv[3] - fwd[3]) < 0.001); + // The Finnish town of Tornio is inside the "another_test" grid + let tio = Coor4D::geo(65.85, 24.10, 0., 0.); + let tio = ellps.cartesian(&tio); + let mut data = [tio]; + ctx.apply(op, Fwd, &mut data)?; + let fwd = data[0]; + dbg!(fwd); + assert!(fwd[0].is_finite()); + + // The Norwegian town of Longyearbyen is outside of both grids + let lyb = Coor4D::geo(78.25, 15.5, 0., 0.); + let lyb = ellps.cartesian(&lyb); + let mut data = [lyb]; + ctx.apply(op, Fwd, &mut data)?; + assert!(data[0][0].is_nan()); + Ok(()) } } diff --git a/src/inner_op/gridshift.rs b/src/inner_op/gridshift.rs index 0bd31be..1168759 100644 --- a/src/inner_op/gridshift.rs +++ b/src/inner_op/gridshift.rs @@ -128,16 +128,25 @@ pub fn new(parameters: &RawParameters, ctx: &dyn Context) -> Result { let def = ¶meters.definition; let mut params = ParsedParameters::new(parameters, &GAMUT)?; - for grid_name in params.texts("grids")?.clone() { - if grid_name.ends_with("@null") { - params.boolean.insert("null_grid"); - continue; + for mut grid_name in params.texts("grids")?.clone() { + let optional = grid_name.starts_with('@'); + if optional { + grid_name = grid_name.trim_start_matches('@').to_string(); } - // TODO: Handle @optional grids + if grid_name == "null" { + params.boolean.insert("null_grid"); + break; // ignore any additional grids after a null grid + } - let grid = ctx.get_grid(&grid_name)?; - params.grids.push(grid); + match ctx.get_grid(&grid_name) { + Ok(grid) => params.grids.push(grid), + Err(e) => { + if !optional { + return Err(e); + } + } + } } let fwd = InnerOp(fwd); @@ -237,6 +246,46 @@ mod tests { Ok(()) } + + #[test] + fn optional_grid() -> Result<(), Error> { + let mut ctx = Plain::default(); + let op = ctx.op("gridshift grids=@../../geodesy/datum/test_subset.datum, @missing.gsb, ../../geodesy/datum/test.datum")?; + + // Copenhagen is outside of the (optional, but present, subset grid) + let cph = Coor4D::geo(55., 12., 0., 0.); + let mut data = [cph]; + + ctx.apply(op, Fwd, &mut data)?; + let res = data[0].to_geo(); + assert!((res[0] - 55.015278).abs() < 1e-6); + assert!((res[1] - 12.003333).abs() < 1e-6); + + ctx.apply(op, Inv, &mut data)?; + assert!((data[0][0] - cph[0]).abs() < 1e-10); + assert!((data[0][1] - cph[1]).abs() < 1e-10); + + // Havnebyen (a small town with a large geodetic installation) is inside the subset grid + let haby = Coor4D::geo(55.97, 11.33, 0., 0.); + let mut data = [haby]; + let expected_correction = Coor4D([11.331, 55.971, 0., 0.]); + ctx.apply(op, Fwd, &mut data)?; + let correction = ((data[0] - haby) * Coor4D([3600., 3600., 3600., 3600.])).to_degrees(); + dbg!(correction); + assert!((correction - expected_correction)[0].abs() < 1e-6); + assert!((correction - expected_correction)[1].abs() < 1e-6); + + Ok(()) + } + + #[test] + fn missing_grid() -> Result<(), Error> { + let mut ctx = Plain::default(); + let op = ctx.op("gridshift grids=missing.gsb"); + assert!(op.is_err()); + + Ok(()) + } } // See additional tests in src/grid/mod.rs