[PROJ] Proposal for a geodetic deformation model format

[Lesparre, Jochem]

Hi,



Kristian wrote:

> Regarding the iterative inverse deformation operation, it works but it isn't particularly good to be honest.

[...]

> Improving this could be a fun challenge for  some one more mathematically inclined than me.





I would like to add two theoretical issues:



1. There is no guarantee that the iterative inverse grid shift operation will always convert to a solution.



The obvious case is when the forward grid shift operation shifts the coordinates of the point of interest outside the bounds of the grid. This often happens near two of the four grid edges. In this case, the inverse operation can't find a solution. This case is handled by PROJ. But I'm not sure about more subtle cases.



During the iteration, the coordinates of the point of interest get shifted multiple times. If in this process the coordinates are shifted back to a previous location, the process will get in an infinite loop. In this loop a point is endlessly shifted back and forth, or in a more complicated pattern like a circle. Of course, such behaviour could not occur in a sensible grid. But most creators of a grid shift file will not check their grid on this.



Another, more likely reason for the inverse grid shift operation not to be able to find a solution is when during the iteration process the coordinates of the point of interest are shifted (temporarily) outside the bounds of the grid.



It could be useful to be able to check grids on these cases?





2. In a NTv2 grid shift file with a subgrid, discontinuities could occur at the edge between the parent grid and the subgrid.



I will explain this with an example. Take a look at the small grid below, with a parent grid (grid points marked with x) and a subgrid (additional grid points marked with +). To interpolate the grid shift at a point of interest on the edge (e.g. point marked with o), one could do this based on the four values of the parent grid points (x) north of the point of interest, or based on the four values of the subgrid points (+) south of the point of interest. Often, the grid shift interpolated from the parent grid and the grid shift interpolated subgrid for a point at the edge (o) will be different. This could cause a discontinuity, as a point slightly north of the point of interest gets a different grid shift as a point slightly south of the point of interest. For a sensible grid, the difference will be so small that it is not significant, but there are no guarantees. A way to check a grid on this could be useful?



It is also possible to prevent such discontinuities to occur at all: by first interpolating two auxiliary points (marked with . ) and then interpolating points with the these two auxiliary points (.) and the points of the subgrid (+). Maybe PROJ. already does interpolate in that way? I didn’t check the source code.



x     x     x     x





x . . x     x     x





x +o+ x + + x     x

+ + + + + + +

+ + + + + + +

x + + x + + x     x



Regards, Jochem



PS: This theoretical issue at the grid edge is not related to the practical issue a grid edges that is reported here: https://github.com/OSGeo/PROJ/issues/1663





-----Original Message-----
From: PROJ <proj-bounces at lists.osgeo.org> On Behalf Of Kristian Evers
Sent: vrijdag 6 december 2019 14:42
To: Even Rouault <even.rouault at spatialys.com>; Chris Crook <ccrook at linz.govt.nz>
Cc: 'proj at lists.osgeo.org' <proj at lists.osgeo.org>
Subject: Re: [PROJ] Proposal for a geodetic deformation model format



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

_______________________________________________

PROJ mailing list

PROJ at lists.osgeo.org<mailto:PROJ at lists.osgeo.org>

https://lists.osgeo.org/mailman/listinfo/proj


Disclaimer:
De inhoud van dit bericht is uitsluitend bestemd voor geadresseerde.
Gebruik van de inhoud van dit bericht door anderen zonder toestemming van het Kadaster 
is onrechtmatig. Mocht dit bericht ten onrechte bij u terecht komen, dan verzoeken wij u 
dit direct te melden aan de verzender en het bericht te vernietigen. 
Aan de inhoud van dit bericht kunnen geen rechten worden ontleend.

Disclaimer:
The content of this message is meant to be received by the addressee only.
Use of the content of this message by anyone other than the addressee without the consent 
of the Kadaster is unlawful. If you have received this message, but are not the addressee, 
please contact the sender immediately and destroy the message.
No rights can be derived from the content of this message
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20191208/5a020886/attachment-0001.html>