diff --git a/.gitignore b/.gitignore index b1749f60..ee8c4487 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,4 @@ /.vscode Cargo.lock *.code-workspace +draft_commit_message.md diff --git a/CHANGELOG b/CHANGELOG index f7801daf..4439f537 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -1,3 +1,128 @@ +2023-11-20: Update code examples and Cargo.toml to reflect v0.11.0 (Thomas Knudsen) +2023-11-20: Remove an unintended version id in the Docs URL (Thomas Knudsen) +2023-11-20: Plain Context improvements (Thomas Knudsen) +2023-11-20: Update kp documentation (Thomas Knudsen) +2023-11-20: Documentation for the deformation operator (Thomas Knudsen) +2023-11-20: a few extra steps in the HOWTO-RELEASE.md recipe (Thomas Knudsen) +2023-11-20: Remove superfluous NTv2 README (Thomas Knudsen) +2023-11-20: .gitignore draft commit message (Thomas Knudsen) +2023-11-17: use once_cell to avoid unsafe blocks (Sean Rennie) +2023-11-16: Grid provider for Plain (#77) (Thomas Knudsen) +2023-11-14: Clean up lint in tests (#76) (Thomas Knudsen) +2023-11-14: Add NTv2 multi subgrid support (#74) (Sean Rennie) +2023-11-10: Merge pull request #75 from busstoptaktik/gridshift (Thomas Knudsen) +2023-11-10: Add new NTv2 grids for testing (Thomas Knudsen) +2023-11-10: Repair grid interpolation problem (Thomas Knudsen) +2023-11-10: works (Thomas Knudsen) +2023-11-10: gridshifting (Thomas Knudsen) +2023-11-10: remove feature gates (#73) (Sean Rennie) +2023-11-09: A number of improvements to the grid handling (#72) (Thomas Knudsen) +2023-11-09: Get NTv2 working + wrap up (Thomas Knudsen) +2023-11-09: Clean up after rebase (Thomas Knudsen) +2023-11-03: add test fixtures (Sean Rennie) +2023-11-03: rework ntv2grid to work with multi grid support (Sean Rennie) +2023-11-03: use within when checking a coord is contained in a grid (Sean Rennie) +2023-11-03: merge fixes (Sean Rennie) +2023-10-20: feature gate ntv2 (Sean Rennie) +2023-10-20: update to work with multi grid test (Sean Rennie) +2023-09-18: somewhat simplify the interpolation code (Sean Rennie) +2023-08-30: add ntv2 grid parser and struct (Sean Rennie) +2023-10-20: Using vector approach to grids (Sean Rennie) +2023-10-17: use multiple grids in gridshift and deformation (Sean Rennie) +2023-10-17: managed getting a GridTrait grid (Sean Rennie) +2023-11-08: Rg proj comparison (#71) (Thomas Knudsen) +2023-11-08: Merge pull request #70 from busstoptaktik/better-basegrid (Thomas Knudsen) +2023-11-08: BaseGrid handles all grid orientations (Thomas Knudsen) +2023-11-07: Better names for Grid trait methods (Thomas Knudsen) +2023-11-03: Remove superfluous assignment in omerc (Thomas Knudsen) +2023-11-03: Minor syntax cleanup in CI (Thomas Knudsen) +2023-11-03: Clippy and fmt in CI + resolve clippy warnings (Thomas Knudsen) +2023-11-03: 67 handle optional grids (#68) (Sean Rennie) +2023-11-01: Document the meaning of `within` (Thomas Knudsen) +2023-11-01: Prototype multi grid support (#65) (Sean Rennie) +2023-10-25: standardise formatting comment (#66) (Sean Rennie) +2023-10-21: Implement Swiss Oblique Mercator (somerc) operator (#63) (Sean Rennie) +2023-10-21: Formatting (Thomas Knudsen) +2023-10-20: Add partial support PROJ ellipsoid definitions (#64) (Sean Rennie) +2023-08-21: Grammar (Thomas Knudsen) +2023-08-21: Avoid double correction for lat_0 in tmerc (Thomas Knudsen) +2023-08-17: Additional Rumination #3 clean up (Thomas Knudsen) +2023-08-17: Clean up of Rumination #3 (Thomas Knudsen) +2023-08-16: Improve documentation by extension, editing, and code restructuring (Thomas Knudsen) +2023-08-15: Brush up error messages (Thomas Knudsen) +2023-08-14: Let proj_parse reject obviously invalid pipeline invocations (Thomas Knudsen) +2023-08-01: Linguistics (Thomas Knudsen) +2023-08-14: Refactor math and preludes, mostly to improve the doc structure (Thomas Knudsen) +2023-08-14: Restructure math and authoring modules (Thomas Knudsen) +2023-08-11: Simplify ParsedParameters (#58) (Thomas Knudsen) +2023-08-10: kp: inv + roundtrip OK (Thomas Knudsen) +2023-08-08: Text series option for ParsedParameters (Thomas Knudsen) +2023-08-07: Jacobian: Remove unused factors (Thomas Knudsen) +2023-08-07: Jacobian: Minor correction (Thomas Knudsen) +2023-08-07: Handle pipeline inversions and omissions (Thomas Knudsen) +2023-08-07: parse_proj: handle a few more pathological cases (Thomas Knudsen) +2023-08-06: Merge pull request #57 from busstoptaktik/proj-parser (Thomas Knudsen) +2023-08-06: PROJ-to-Rust Geodesy syntax converter (Thomas Knudsen) +2023-08-04: Jacobian (#54) (Thomas Knudsen) +2023-08-01: Restore support for inline comments in step definitions (Thomas Knudsen) +2023-08-01: Merge pull request #53 from busstoptaktik/assertions (Thomas Knudsen) +2023-08-01: Better assertions (Thomas Knudsen) +2023-08-01: Merge pull request #52 from busstoptaktik/introspection (Thomas Knudsen) +2023-08-01: Operator introspection functionality (Thomas Knudsen) +2023-08-01: Spelling (Thomas Knudsen) +2023-08-01: Streamline split_into_steps() (Thomas Knudsen) +2023-07-30: Better exposure of Coor2D and Coor3D (Thomas Knudsen) +2023-07-30: Move proj operator to examples (Thomas Knudsen) +2023-07-11: Linting: cargo clippy + cargo fmt (Thomas Knudsen) +2023-07-11: Correction of an example (Thomas Knudsen) +2023-07-11: Remove the unused 'assets' tree (Thomas Knudsen) +2023-07-11: Remove file accidentally committed in #6e58f43 (Thomas Knudsen) +2023-07-11: Rumination #2: Linguistics and mino corrections (Thomas Knudsen) +2023-07-09: Rename the NMEA operators to dm/dms and iso_dm/iso_dms (Thomas Knudsen) +2023-07-09: Rename support files from 'collection' to 'register' (Thomas Knudsen) +2023-07-09: Experiment with a different formatting convention (Thomas Knudsen) +2023-07-09: Add missing half sentence to nkg.collection (Thomas Knudsen) +2023-07-09: Extend the 'plain' context provider and use it! (Thomas Knudsen) +2023-07-08: Substantial improvements of the 'plain' context provider (Thomas Knudsen) +2023-07-08: Merge pull request #44 from busstoptaktik/kp (Thomas Knudsen) +2023-07-08: Merge branch 'main' into kp (Thomas Knudsen) +2023-07-08: Remove some superfluous messages from tmerc (Thomas Knudsen) +2023-07-08: Switch examples from simple-logger to env-logger (Thomas Knudsen) +2023-07-08: parse_sexagesimal now supports NaN (Thomas Knudsen) +2023-07-08: KP overhaul (Thomas Knudsen) +2023-07-08: Reintroduce 'clap' as de-dependency (Thomas Knudsen) +2023-07-07: Adjustment of dependencies (Thomas Knudsen) +2023-07-06: Merge pull request #43 from busstoptaktik/tmerc (Thomas Knudsen) +2023-07-06: Less verbosity for tmerc (Thomas Knudsen) +2023-07-06: Rough overhaul (Thomas Knudsen) +2023-07-05: Remove intentional test failure (Thomas Knudsen) +2023-07-04: Improved description of deformation based transformations (Thomas Knudsen) +2023-07-02: Forward part of the deformation operator (Thomas Knudsen) +2023-07-02: Syntax correction in example (Thomas Knudsen) +2023-07-01: Improved code example in rumination 000 (Thomas Knudsen) +2023-07-01: Merge pull request #42 from busstoptaktik/coord_types (Thomas Knudsen) +2023-07-01: Correct a formal error in a comment (Thomas Knudsen) +2023-07-01: Implement CoordinateSet for tuples (T, f64, f64) & (T, f64) (Thomas Knudsen) +2023-06-30: Add remaining examples to justfile run-all (Thomas Knudsen) +2023-06-27: Minor cleanups (Thomas Knudsen) +2023-06-27: Enter Coor3D (Thomas Knudsen) +2023-06-27: Converters for Coor32 (Thomas Knudsen) +2023-06-27: More coordinate mopping (Thomas Knudsen) +2023-06-27: Rename coord.rs to coor4d.rs (Thomas Knudsen) +2023-06-27: Mop up after rename (Thomas Knudsen) +2023-06-26: Rename Coord to Coor4D (Thomas Knudsen) +2023-06-26: Minor improvements to the slice-example (Thomas Knudsen) +2023-06-26: Give an example of using slices [Coord] (Thomas Knudsen) +2023-06-24: Clean up (Thomas Knudsen) +2023-06-24: Correction of some remaining old parameter syntax (Thomas Knudsen) +2023-06-23: Example of user defined coordinate types and containers (Thomas Knudsen) +2023-06-23: Exit DirectPosition (Thomas Knudsen) +2023-06-23: Eliminate the need for Plain context in tests (Thomas Knudsen) +2023-06-07: Ironing out some 'build with_plain' glitches (Thomas Knudsen) +2023-06-07: Minor addtion to HOWTO-RELEASE (Thomas Knudsen) +2023-06-07: Split dependencies between lib and cli (Thomas Knudsen) +2023-06-07: HOWTO-RELEASE. Draft but functionally complete (Thomas Knudsen) +2023-06-07: More descriptive keywords (Thomas Knudsen) 2023-06-07: First preparations for 0.10.0 (Thomas Knudsen) 2023-06-07: Linguistics (Thomas Knudsen) 2023-06-06: UTM: Support "south" aspect (Thomas Knudsen) diff --git a/CHANGES.md b/CHANGES.md new file mode 100644 index 00000000..575aa5d9 --- /dev/null +++ b/CHANGES.md @@ -0,0 +1,36 @@ +# Geodesy + +## 0.11.0 Release notes + +### Improvements + +- Handle lists-of-grids, `@optional` grids, and the `@null` grid in `grids=` clauses +- Support NTv2 format datum shift grids +- Overall documentation brush up and extension +- Implement `somerc`, the Swiss Oblique Mercator operator +- Implement `deformation`, the 3D intrapalte deformation operator +- Rename the `NMEA` operators to `dm/dms` and `iso_dm/iso_dms` +- Support jacobian-of-projection and the corresponding deformation factors +- `proj_parse`: Translate PROJ syntax to Rust Geodesy syntax + with partial support for PROJ ellipsoid definitions +- Through `proj_parse`, the `Plain` context, and hence `kp` now supports PROJ syntax + (although with the limitations implied by Geodesy only supporting parts of the PROJ + operator gamut, and by the differing input-output conventions between `kp` and `proj/cct`) + as demonstrated by this example: + + ```sh + echo 55 12 | kp geo:in | kp "proj=pipeline step proj=utm zone=32 step inv proj=utm zone=32" | kp geo:out + ``` + + which does nothing, in a very convoluted way +- Partial operator introspection +- General support for 2D, 3D, 4D, and 32 bit 2D coordinates (`Coor2D`, `Coor3D`, `Coor4D`, `Coor32`) +- Hence `Coord` is gone + +### Bug fixes + +- Avoid double correction for lat_0 in inverse tmerc + +### Acknowledgements + +A huge thank you goes to [Sean Rennie](https://github.com/Rennzie) who a.o. did most of the work on the improved grid support in 0.11.0, and to [Kyle Barron](https://github.com/kylebarron) for pushing Geodesy over the WASM barrier in 0.10.0 diff --git a/Cargo.toml b/Cargo.toml index af6d93d5..4ac677fd 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -4,7 +4,7 @@ description = "A platform for experiments with geodetic transformations and data keywords = ["geospatial", "geodesy", "cartography", "geography"] categories = ["science"] license = "MIT OR Apache-2.0" -version = "0.10.0" +version = "0.11.0" rust-version = "1.65" # Uses the let-else construct introduce in 1.65 authors = ["Thomas Knudsen "] readme = "README.md" diff --git a/HOWTO-RELEASE.md b/HOWTO-RELEASE.md index 5e49a86e..612b0d16 100644 --- a/HOWTO-RELEASE.md +++ b/HOWTO-RELEASE.md @@ -1,9 +1,15 @@ # DRAFT: How to release -- update `Cargo.toml` with new version id, e.g. `"0.11.0"` +- manually check that all [issues](https://github.com/busstoptaktik/geodesy/issues/) + assigned to the + [milestone for upcomming release](https://github.com/busstoptaktik/geodesy/issues?q=is%3Aopen+is%3Aissue+milestone%3A0.11.0) + are resolved +- update `Cargo.toml` and `README.md` with new version id, i.e. `"0.11.0"` + - `just check-all` -- `just changes` preview a new `CHANGELOG` -- `just changelog` generate a new `CHANGELOG` +- `just changes` (to preview a new `CHANGELOG`) +- update `CHANGES.md` +- `just changelog` (to generate a new `CHANGELOG`) - `git commit ...` - `git push` - `git tag v0.11.0` @@ -13,3 +19,6 @@ - `git push ...` - `git switch main` - `cargo publish` +- update `HOWTO-RELEASE.md` to say 0.12 +- `git commit -a -m "Start of work towards 0.12.0` +- `git push ...` diff --git a/README.md b/README.md index 68030541..06179780 100644 --- a/README.md +++ b/README.md @@ -25,7 +25,7 @@ Add Geodesy to your `Cargo.toml` file ```toml [dependencies] -geodesy = "0.10.0" +geodesy = "0.11.0" ``` Then try this minimal example, computing the UTM coordinates of some Scandinavian capitals diff --git a/ruminations/002-rumination.md b/ruminations/002-rumination.md index 3bbc2285..8a93b613 100644 --- a/ruminations/002-rumination.md +++ b/ruminations/002-rumination.md @@ -6,7 +6,7 @@ Thomas Knudsen Sean Rennie -2021-08-20. Last [revision](#document-history) 2023-11-02 +2021-08-20. Last [revision](#document-history) 2023-11-20 ### Abstract @@ -24,6 +24,8 @@ $ echo 553036. -124509 | kp "dms:in | geo:out" - [`adapt`](#operator-adapt): The order-and-unit adaptor - [`cart`](#operator-cart): The geographical-to-cartesian converter - [`curvature`](#operator-curvature): Radii of curvature +- [`deformation`](#operator-deformation): Kinematic datum shift using a + 3D deformation model in ENU-space - [`dm`](#operator-dm): DDMM.mmm encoding. - [`dms`](#operator-dms): DDMMSS.sss encoding. - [`geodesic`](#operator-geodesic): Origin, Distance, Azimuth, Destination and v.v. @@ -186,6 +188,146 @@ curvature prime ellps=GRS80 --- +### Operator `deformation` + +**Purpose:** +Kinematic datum shift using a 3D deformation model in ENU-space + +**Description:** + +Based on Kristian Evers' implementation of the +[corresponding PROJ operator](https://github.com/OSGeo/PROJ/blob/effac63ae5360e737790defa5bdc3d070d19a49b/src/transformations/deformation.cpp). +The deformation operation takes cartesian coordinates as input and +yields cartesian coordinates as output. The deformation model is +assumed to come from a 3 channel grid of deformation velocities, +with the grid georeference given as geographical coordinates in a +compatible frame. + +#### The Deformation + +The deformation expressed by the grid is given in the local +east-north-up (ENU) frame. It is converted to the cartesian XYZ +frame when applied to the input coordinates. +The total deformation at the position P: (X, Y, Z), at the time T1 is +given by: + +```txt + DX(X, Y, Z) = (T1 - T0) * Vx(φ, λ) + (1) DY(X, Y, Z) = (T1 - T0) * Vy(φ, λ) + DZ(X, Y, Z) = (T1 - T0) * Vz(φ, λ) +``` + +where: + +- (X, Y, Z) is the cartesian coordinate tuple for P +- (DX, DY, DZ) is the deformation along the cartesian earth centered + axes of the input frame +- (Vx, Vy, Vz) is the deformation velocity vector (m/year), obtained + from interpolation in the model grid, and converted from the local + ENU frame, to the global, cartesian XYZ frame +- (φ, λ) is the latitude and longitude, i.e. the grid coordinates, + of P, computed from its cartesian coordinates (X, Y, Z) +- T0 is the frame epoch of the kinematic reference frame associated + with the deformation model. +- T1 is the observation epoch of the input coordinate tuple (X, Y, Z) + +#### The transformation + +While you may obtain the deformation vector and its Euclidean norm +by specifying the `raw` option, that is not the primary use case for +the `deformation` operator. Rather, the primary use case is to *apply* +the deformation to the input coordinates and return the deformed +coordinates. Naively, but incorrectly, we may write this as + +```txt + X' = X + DX = X + (T1 - T0) * Vx(φ, λ) + (2) Y' = Y + DY = Y + (T1 - T0) * Vy(φ, λ) + Z' = Z + DZ = Z + (T1 - T0) * Vz(φ, λ) +``` + +Where (X, Y, Z) is the *observed* coordinate tuple, and (X', Y', Z') +is the same tuple after applying the deformation. While formally +correct, this is not the operation we intend to carry out. Neither +are the names used for the two types of coordinates fully useful +for understanding what goes on. + +Rather, when we transform a set of observations, we want to obtain the +position of P at the time T0, i.e. at the *epoch* of the deforming +frame. In other words, we want to remove the deformation effect such +that *no matter when* we go and re-survey a given point, we will always +obtain the same coordinate tuple, after transforming the observed +coordinates back in time to the frame epoch. Hence, for the forward +transformation we must *remove* the effect of the deformation by negating +the sign of the deformation terms in eq. 2: + +```txt + X' = X - DX = X - (T1 - T0) * Vx(φ, λ) + (3) Y' = Y - DY = Y - (T1 - T0) * Vy(φ, λ) + Z' = Z - DZ = Z - (T1 - T0) * Vz(φ, λ) +``` + +In order to be able to discuss the remaining intricacies of the task, we +now introduce the designations *observed coordinates* for (X, Y, Z), and +*canonical coordinates* for (X', Y', Z'). + +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* +(this is necessary, since the deforming system changes as time goes). + +But we cannot *observe* anything with respect to the canonical system: +It represents the world as it was at the epoch of the system. So the observed +coordinates are given in a system slightly different from the canonical. +The deformation model makes it possible to *predict* the coordinates we will +observe at any given time, for any given point that was originally observed +at the epoch of the system. + +But we are really more interested in the opposite: To look back in time and +figure out "what were the coordinates at time T0, of the point P, which we +*actually observed at time T1*". + +But since the georefererence of the deformation grid is given in the canonical +system, we actually need to know the canonical coordinates already in order to +look up the deformation needed to convert the observed coordinates to the +canonical, leaving us with a circular dependency ("to understand recursion, we +must first understand recursion"). + +To solve this, we do not actually need recursion - there is a perfectly +fine solution based on iteration, which is widely used in the inverse case +of plain 2D grid based datum shifts (whereas here, we need it in the forward +case). + +There is however an even simpler solution to the problem - simply to ignore it. + +The deformations are typically so small compared to the grid node distance, +that the iterative correction is way below the accuracy of the transformation +grid information, so we may simply look up in the grid using the observed +coordinates, and correct the same coordinates with the correction obtained +from the grid. + +For now, this is the solution implemented here. + +| Parameter | Description | +|-----------|-------------| +| `inv` | Inverse operation: output-to-input datum. Currently implemented using sign reversion, *without* iterative refinement | +| `raw` | Replace the input coordinate by the correction values, rather than applying them | +| `dt` | Specify a fixed deformation interval, rather than using the difference between `t_epoch` and the point coordinate time | +| `t_epoch` | The temporal origin of the deformation proces, given as decimal year | +| `ellps` | The ellipsoid for the deforming system. Used for converting the ENU elements of the grid, to dLat, dLon, dHeight corrections | +| `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 hence do block instantiation of the operator if they are unavailable. Additionally, if the `@null` parameter is specified as the last grid, points outside of the grid coverage will be passed through unchanged, rather than being stomped on with the NaN shoes and counted as errors | + +**Example**: + +```txt +deformation dt=1000 ellps=GRS80 grids=test.deformation + +deformation raw dt=1000 grids=test.deformation,@another.deformation,@null +``` + +**See also:** The documentation for the corresponding [PROJ operator](https://proj.org/en/9.3/operations/transformations/deformation.html) + +--- + ### Operator `dm` **Purpose:** Convert from/to the ISO-6709 DDDMM.mmm format. @@ -305,7 +447,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 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. | +| `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 hence do block instantiation of the operator if they are unavailable. Additionally, if the `@null` parameter is specified as the last grid, points outside of the grid coverage will be passed through unchanged, rather than being stomped on with the NaN shoes and counted as errors | The `gridshift` operator has built in support for the **Gravsoft** grid format. Support for additional file formats depends on the `Context` in use. @@ -745,3 +887,4 @@ Major revisions and additions: - 2023-07-09: dm and dms liberated from their NMEA overlord - 2023-10-19: Add `somerc` operator description - 2023-11-02: Update `gridshift` operator description with multi, optional and null grid support +- 2023-11-20: Add documentation for the `deformation` operator diff --git a/ruminations/003-rumination.md b/ruminations/003-rumination.md index b2f629e5..adde94d8 100644 --- a/ruminations/003-rumination.md +++ b/ruminations/003-rumination.md @@ -4,7 +4,7 @@ Thomas Knudsen -2021-08-28. Last [revision](#document-history) 2022-05-08 +2021-08-28. Last [revision](#document-history) 2023-11-20 ### Abstract @@ -88,39 +88,28 @@ The `help` option gives the list of options: ```txt $ kp --help -kp 0.7.1 -KP: The Rust Geodesy "Coordinate Processing" program. -Called `kp` in honor of Knud Poder (1925-2019), the -nestor of computational geodesy, who would have found -it amusing to know that he provides a reasonable -abbreviation for something that would otherwise have -collided with the name of the Unix file copying program `cp` - -USAGE: - kp.exe [FLAGS] [OPTIONS] [FILE]... - -FLAGS: - -d, --debug Activate debug mode - -e, --echo Echo input to output - -h, --help Prints help information - -i, --inv Inverse. Use of `inverse` mode excludes the use - of `roundtrip` mode - -r, --roundtrip Roundtrip mode - a signature feature of Knud - Poder's programs: Evaluate the accuracy of the - transformation by comparing the input argument - with its supposedly identical alter ego after - a forward+inverse transformation pair. Use of - `roundtrip` mode excludes the use of `inverse` - mode - -V, --version Prints version information - -v, --verbose Verbose mode (-v, -vv, -vvv, etc.) - -OPTIONS: - -o, --output Output file, stdout if not present - -ARGS: - Operation to apply - ... Files to process + +KP: The Rust Geodesy 'Coordinate Processing' program + +Usage: kp.exe [OPTIONS] [ARGS]... + +Arguments: + The operation to carry out e.g. 'kp "utm zone=32"' + [ARGS]... The files to operate on + +Options: + --inv Inverse operation + -z, --height Specify a fixed height for all coordinates + -t, --time