[PROJ] Proposal for a geodetic deformation model format

[Brian Shaw]

Thanks Evan for the questions and Chris for the answers.

I know that at NGS we are working to solve how to implement many of the 
same challenges as LINZ except for all the US and its territories.  NGS 
plans to essentially have two models to handle deformation with an Inter 
Frame Velocity Model (IFVM) handling horizontal shifts and a Dynamic 
Geoid (DGEOID) to handle geopotential (vertical) shifts.  I know that 
New Zealand deals with earthquakes a bit more than we do but California, 
Alaska and the Pacific have many interesting events and challenges.  
Earlier this year the China Lake earthquake in California moved the 
local surface 12 feet or more in some places but thankfully it was not 
in a populated area with significant infrastructure.  Such major event 
shifts occasionally happen mostly in California, Alaska or the Pacific 
but we are working to create a system that allows us to handle such 
challenges.

Either way Im glad this is getting discussed since in the future system 
of the US time of observations or "epoch" will be essential and the 
coordinates you collect will be good for the day of observation.  The 
new paradigm acknowledges that the world is dynamic and constantly 
changing and so will coordinates over time. This will make the metadata 
of when a coordinate was determined *essential*, particularly if you are 
trying to forecast where a physical place is moving or where it has been.

While I havent got my colleagues to embrace the term since its a 
mouthful, I have been calling the future coordinates spatiotemporal 
coordinates.  I think we at NGS need to talk a lot more with LINZ about 
the grids, models and even database schemas.

Cheers
Brian

On 12/5/2019 5:50 PM, Chris Crook wrote:
> Hi Even
>
> Thanks for detailed reading and questions.
>
>> - regarding "The horizontal displacement is represented by metres east and
>> north. Typically base datum will be a geographic coordinate system using
>> latitude and longitude and it would be more efficient to define displacement
>> in terms of degrees.".
>> Is it the same convention as PROJ deformation method
>> (https://proj.org/operations/transformations/deformation.html) which uses
>> grids to store corrections in the ENU space as far as I understand ?
>>
>> (more a question for Kristian: I presume this ENU space is the same as the
>> UVW topocentric coordinate system described in 4.1.2
>> "Geocentric/topocentric conversions" of https://www.iogp.org/wp-
>> content/uploads/2019/09/373-07-02.pdf ?)
>>
>> For the easting, northing I presume yes. But I'm not sure if the convention
>> you mention for vertical correction "The vertical direction will be the local
>> vertical of the coordinate system." corresponds to the U of ENU. Or perhaps
>> this falls in your remark "For deformation there is no practical need to
>> distinguish gravitational and geometric vertical" ?
>> That said, PROJ deformation method is a bit different in that it is meant to
>> perate with input and output coordinates in geocentric space: it uses
>> internally a geocentric to geodetic conversion to be able to retrieve the
>> correction values from the grid(s), and then convert those corrections from
>> the ENU space to the geocentric space.
> I hadn't realised PROJ would be working in a geocentric space.
>
> What I was thinking when I suggested degrees would be more efficient was that
> had converted coordinates to geographic lon/lat/ellipsoidal height (as you need lat/lon
> to interpolate), then applying a correction to the lon/lat coordiantes in metres requires (a little)
> more calculation than applying corrections in degrees.   This could make more difference
> for an inverse calculation, as if you are iterating you could end up converting multiple times
> between geocentric and geographic (and geocentric to geographic is itself iterative).
>
> However if we are applying the correction to geocentric coordinates then
> this doesn't apply.
>
> For the inverse calculation it is probably sensible to work in the geographic coordinates
> if the solution is iterated, and then convert to geocentric once the inverse
> deformation has been calculated.
>
>> - regarding "Strictly the inverse calculation of the deformation model is an
>> iterative calculation, as the location of a point in the interpolation coordinate
>> system is not known until it is transformed back to that system."
>> I just want to mention that PROJ uses an iterative method for the inverse of
>> horizontal shift grids (NTv2 etc), so wouldn't be a problem to use that for
>> deformation models as well. Actually PROJ deformation method also uses
>> that iterative method.
> The NTv2 grids were originally designed for transformations between quite
> different datums that pre-dated satellite geodesy, and for which coordinate
> changes of a few hundred metres are quite likely, so in that case the iteration
> could be quite important.
>
> The coordinate changes due to deformation are likely to be much smaller (hopefully!).
> I think 10 metres would be an upper bound.  On the other hand the deformation will
> vary more rapidly with location in the event of earthquake related deformation.  I'll
> do some testing with the LINZ deformation model to see what difference it makes.
>
> But good to know that PROJ will do it the right way :-)
>
>> - regarding nesting of grids, is my understanding correct that if you have a
>> grid G, that has subgrids G_1 et G_2, and that a point falls in G and G_1, then
>> the corections of both G and G_1 should be applied ? Just wanted to
>> underline this, as it is different of the NTv2 model where you apply one or
>> the other (preferably G_1), but not both.
> No - I've tried to make this clearer in the document.  The proposal is to use the
> same approach as eg NTv2.   A single component (eg deformation due to an
> earthquake) may be defined by multiple grids.  At any location only one of those
> grids will apply.
>
> However the deformation model as a whole has multiple components (eg
> secular velocity, earthquake 1, earthquake 2, ...) and these are evaluated
> separately and summed together.
>
>> - it could help for a better understanding to have a small example of a CSV
>> file illustrating the different concepts you expose. Nested grids, different
>> time functions etc.
> Good idea - I'll build a CSV file matching the NZGD2000 deformation model.
>
>> - The approach master CSV file + separate grids in dedicated files has its
>> merits indeed. The actual grid format, GTiff, HDF5 or whatever raw binary
>> format that can support multiple samples per grid node will make little
>> difference because the specificities/complexities of the description of the
>> deformation model is in the CSV file itself.
>>
>> - So if a earthquake happens, basically you'll create a new grid in the area that
>> is affected, and create a new version of the CSV file where you add an extra
>> row that references it, right ? (I'm probably over-simplifying)
> That would be the most likely result.  For the last earthquake Kaikoura there
> would be two rows, as we treated near field and far field deformation
> differently (we still need to write this up).  Also we had two releases of the
> deformation model.  The second release was about a year later once we
> had better information about the total deformation (the update is an awkward
> and largely inseparable mix of better estimation of co-seismic deformation
> and additional post-seismic deformation).
>
>> - Any indication of the number of grids a deformation model might have ? For
>> the NZ case for example
> The NZGD2000 has a secular velocity grid plus grids for 12 events since 2000.0
> 124 grids at last  count I think, but that is individual grids, not the number of nested grid.
> Each event could be represented by 2 or 3 components with different time
> functions, so probably between 20 and 40 in total.
>
> However as our ability to measure and model deformation quickly is improving
> then the rate of  generating grids increases.
>
>> - Just a detail: regarding "velocity f(t) = (t – t0)/(365.2425 days)", if time and
>> volicities are expressed respectively in decimal year and
>> some_linear_unit/decimal_year, that should probably be just f(t) = t-t0
> Yes - I copied this directly from the LINZ deformation model documentation.
> Practically this is no different to expressing the time in decimal years.  I haven't
> defined anywhere the way in which date/time are defined, but we use
> a full date/time format rather than decimal years for defining the model as
> that is more intuitive for earthquakes.  And of course most data we are transforming
> has a specific date/time.   So time differences are determined in terms of "real"
> time units such as seconds, days, or whatever.
>
>> - I'm still having issues to grasp the practical implications of the "reverse
>> step". Is my understanding correct that a coordinate computed for example
>> for
>> NZGD2000 @ 2010.0 with a deformation model of 2019 might potentially
>> change with a deformation model of 2020 ? Should the deformation model
>> used be also captured in the metadata associated with a coordinate? If so, it
>> looks like a change in the deformation model would be better handled by
>> defining a revised version of the datum (but I know you don't want to publish
>> a new NZGD datum every 6 month). I'm super confused... I'm pretty sure we
> It absolutely is confusing!   Technically whenever we change the deformation
> model it is a new version of NZGD2000.  However the world, and particular the
> users of NZGD2000 coordinates,  was not ready for this, and possibly still isn't.
> Practically the time required to propagate new coordinate systems through the
> made creating multiple versions of the coordinate system impractical.   Also there
> are 30 projection coordinate systems associated NZGD2000, so we would need
> to create multiple versions of those as well.  The very recent changes to coordinate
> reference system models provide for concepts that will help with this such as
> datum ensembles.  Having PROJ able to retrieve coordinate system definitions
> from EPSG is a huge step towards making this practical in the open source stack.
> Having software that can handle deformation is another huge step.
>
> The way we have been managing NZGD2000 is to try to protect the NZ users
> from dealing with this more than necessary.   This has been very pragmatic.
> Since the tools most users have available are not "deformation aware", we
> have been managing the coordinate system to allow users to ignore deformation
> as far as possible.
>
> As far as possible is to provide a coordinate reference system that hide the
> deformation.  (I personally think NZGD2000 is better described as a geospatial reference
> system than a geodetic datum).
>
> To take the example of the coordinate NZGD2000 at 2010.  When we change
> the deformation model we certainly are changing the relationship between
> an NZGD2000 coordinate and the corresponding ITRF coordinate.
> However we do preserve the relationship  between the NZGD2000 coordinate
> and the physical location that it refers to.  So you don't need to change the
> coordinates stored in your GIS systems etc - they are still correct for the features
> they refer to.
>
> As far as possible doesn't cover very heterogeneous deformation close to
> an earthquake epicentre.  But for example in the case of the 2016 Kaikoura earthquake
> the city of Wellington 150km from teh epicentre was shifted 20cm as a whole.
> This is absorbed in the deformation model so that the NZGD2000 coordinates of marks do not change.
>
> If you have a data set for Wellington observed in 2010 in terms of an ITRF based
> coordinate system and convert to NZGD2000 then it will used the relationship
> between ITRF and NZGD2000 that applied in 2010.  The NZGD2000 coordinate
> that you get is the same as if you had observed the same physical feature in
> terms of ITRF in 2019 it would have a different ITRF coordinate (both because
> of the earthquake and the 5cm/year secular deformation).  However if you
> and convert it to NZGD2000 using the post earthquake version of the deformation
> model you will get the same NZGD2000 coordinate (give or take errors of observation
> and errors in the deformation model).
>
> Nearer the epicentre the situation is different - the deformation is too intense to hide
> from users and so in that area the NZGD2000 coordinates of points do change.  LINZ
> provided grids of the deformation to allow users to update their data sets to account
> for the deformation (as well as updating the deformation model).  The good thing is
> that in this region users expect the coordinates to change, and in fact are asking
> for updates.
>
>> had that discussion already, but I've forgotten the outcome. I'd be very
>> interested in having a concrete example of a workflow of how users are
>> supposed to use such deformation models, and which metadata they are
>> supposed to store to be able to later reprocess their coordinates.
> This is the critical question, and one that we don't have a good answer to!
> At  the moment there is no workflow as virtually no software can use the
> deformation model.  That is why NZGD2000 is implemented the way it is,
> so that NZGD2000 data sets from different epochs can be compared and
> used without having to account for deformation.
>
> The nearest there is to a workflow is the process for users near an area
> affected by deformation where they can use the grids we provide to update
> their data (we provide NTv2 grids for the horizontal component).  We have
> been calling this a "patch".
>
> Your question about metadata is key.  The general answer is that (as for
> all coordinates) they need to be aware of what coordinates are referenced
> to when they are originally captured, and what transformations have been applied to
> them since then.   In areas affected by a patch this is more critical.
>
> Looking to the future this could change in many ways.  We could create
> new versions of the datum for each update of the deformation model.  We
> could imagine that coordinates are held with their effective epoch, and if
> users want to compare data sets in terms of the physical bit of ground they
> refer to they will transform the data sets to a common epoch on the fly.
> Many possibilities!  We are not there yet, but whatever it is the thing we do
> need is the ability to define and apply complex time dependent transformations.
>
> Whew! That was a long answer.  I hope this helps a bit.  At the recent FOSS4G Oceania
> conference I presented a talk on this.   https://www.youtube.com/watch?v=xdvI9Yrs4iA&t=141s
> There is also some information on the LINZ website at https://www.linz.govt.nz/data/geodetic-system/datums-projections-and-heights/geodetic-datums/new-zealand-geodetic-datum-2000-nzgd2000/nzgd2000-deformation-model,
> but most of this is detail rather than conceptual.
>
> One thing we are very aware at LINZ is that we have not had sufficient resource/time
> to publish some of the background to how we have implemented recent deformation
> model updates, but conceptually the driver hasn't changed, which is to keep things
> working as simply as we can for the users of the NZGD2000 coordinates, most of whom are
> probably unaware of the deformation model or the need for one.
>
>> Even
> ________________________________
>
> This message contains information, which may be in confidence and may be subject to legal privilege. If you are not the intended recipient, you must not peruse, use, disseminate, distribute or copy this message. If you have received this message in error, please notify us immediately (Phone 0800 665 463 or info at linz.govt.nz) and destroy the original message. LINZ accepts no responsibility for changes to this email, or for any attachments, after its transmission from LINZ. Thank You.
> _______________________________________________
> PROJ mailing list
> PROJ at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/proj

-- 
*************************************
Brian Shaw
Geodesist
NOAA/NOS/National Geodetic Survey
Phone # 240-533-9522

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20191205/6bced98c/attachment-0001.html>