[PROJ] Proposal for a geodetic deformation model format

[Kristian Evers]

Chris et al,

I am late to the party so let me just try to comment on a few things from the discussion below as well as provide a few use cases where I would benefit from this new grid format. Generally I also think this is a very well thought out and thorough description of a new format. It covers most bases for my needs and I think that the "multi file" solution aligns quite closely with current PROJ practice so it shouldn't be too complicated to make things work. Sorry for the long text - I got carried away, as often happens when the subject is interesting!


COORDINATE SPACE(S)

I originally wrote the deformation operation with the sole purpose of providing a practical implementation of the transformations described in the paper by Häkli et al [0]. This obviously influenced the implementation and there is the possibility that it isn't generic enough for all purposes, although I haven't yet come across a situation that wouldn't fit regarding the use 3D deformation models in coordinate transformations.

The deformation operation in PROJ works in a mix of geocentric and geographic coordinates [1]. Actual adjustments are done in geocentric/cartesian space. This is analogous to e.g. a Helmert operation which also works in geocentric space. And as Even points out the deformation operation is often use in a pipeline in conjunction with one or more Helmert operations. I generally prefer to work in geocentric space when doing this sort of work since it is the "native" coordinate space of the frame and because it is nice to work within a linear coordinate systems with axes that intersects at right angles. Unfortunately it is not really possible to create a grid of surface displacements or velocities in geocentric space without creating a massive sparse 3D matrix which is not practical in any way. This is the reason why the deformation operation requires grids in topocentric/ENU space and internally converts between the two coordinate systems. It does come with a bit of overhead, especially in the inverse mode, but it isn't too bad compared to let's say a 14-parameter Helmert operation.

Doing gridshift adjustments in degrees generally works quite well and it what has been done for years and years in PROJ. The main disadvantage, I think, is that it is difficult for the human eye to interpret values in the grid. I tend to find it quite useful to plot the grids I use in QGIS and having the values expressed in meters or millimeters gives me an immediate advantage. There's nothing in the way to support both types of gridshifts. I do think though that expressing velocities in terms of degrees is a bad idea (I don't think that was the intention, just putting it out there). 

Regarding the iterative inverse deformation operation, it works but it isn't particularly good to be honest. I reworked the existing iterative algorithm used for horizontal grid shifts into a 3D version and I never quite got it right. It doesn't perform well in a round-trip test and after about ten trips back and forth it has introduced too much noise due to numerical inaccuracy in the algorithm. Improving this could be a fun challenge for  some one more mathematically inclined than me.

I agree that it doesn't make much of a difference whether adjustments in the up components are done in relation to a gravimetric or geometric reference.


GRID DELIVERY METHOD (SINGLE VS. MULTI FILE)

I agree with the proposal that the multi file solution is preferable. Especially if zipped up a in aggregate package. This was also touched upon in the original GitHub discussion 18 months ago. I think JSON would be preferable to a CSV based format since it offers a bit more structure and you would be able to define a proper JSON schema that can be used to validate the format of the file. PROJ has the ability to parse JSON already. 

A package of a number of grid files and a metadata file describing the use of the grids should be fairly simple to translate to a PROJ pipeline, provided that the usage of the grids align with operations available in PROJ. It would of course be much simpler if everything was just expressed as a PROJ pipeline but I guess we can't assume that everyone will be using PROJ :-)


RELEVANT USE CASE: GREENLAND

Greenland is subject to significant elastic deformation due to rapid mass-loss from the ice sheet. The deformation varies a lot from year to year which creates a need for a rather complicated deformation model. At the moment I have solved the problem by creating a set deformation models that cover a specific year each. That is, for every year since 1996 I have a 3D deformation model that only applies for that particular year. The plan is to create a pipeline of deformation operations that applies as many adjustments from the deformation models as needed. E.g. if a coordinate from 2010 is transformed it would have to pass through 14 deformation operations (2010-1996) to account for all the deformation. I think for the final transformations I will narrow the number of gridded models down so that they each cover a three year period. The point here is that it would be nice to be able to support this type of transformation in the new deformation model format. If I read the proposal correctly this would be possible but I may have missed some details. PROJ is actually not able to do this as of now but I will extend the deformation operation to allow this in the future (functionality similar to the +t_epoch/+t_final setup in hgridshift is needed). 


RELEVANT USE CASE: ICELAND

Similarly to New Zealand, Iceland exists in very a geodynamically complicated setting. Transformations in Iceland should ideally take care of secular (GIA, divergent plate boundary) and episodic deformation (earth quakes, volcanism) as well as the regular Helmert transformations involving ITRFxxxx. In PROJ this can be handled as a pipeline including the helmert, hgridshift, vgridshift and deformation operations. I gave a talk on the subject last year which covers this use case and some thoughts about a practical implementation [2]. Recently that work described in that presentation has culminated in a pull request for the proj-datumgrid repository [3].


/Kristian


[0] Häkli, P., Lidberg, M., Jivall, L., Nørbech, T., Tangen, O., Weber, M., Pihlak, P., Aleksejenko, I., and Paršeliunas, E. The NKG2008 GPS campaign – final transformation results and a new common Nordic reference frame. Journal of Geodetic Science, 6(1):1–33, 2016. https://doi.org/10.1515/jogs-2016-0001 

[1] Yes, I use formulas equivalent to what is described in section 4.1.2 "Geocentric/topocentric 
conversions" of https://www.iogp.org/wp-content/uploads/2019/09/373-07-02.pdf 

[2] https://www.fig.net/fig2018/rfip/6_RFIPIstanbulKristianEvers.pdf 

[3] https://github.com/OSGeo/proj-datumgrid/pull/67