<html>
  <head>
    <meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
  </head>
  <body>
    Thanks Evan for the questions and Chris for the answers.<br>
    <br>
    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.<br>
    <br>
    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 <b>essential</b>,
    particularly if you are trying to forecast where a physical place is
    moving or where it has been.<br>
    <br>
    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.<br>
    <br>
    Cheers<br>
    Brian<br>
    <br>
    <div class="moz-cite-prefix">On 12/5/2019 5:50 PM, Chris Crook
      wrote:<br>
    </div>
    <blockquote type="cite"
cite="mid:A87E66F06E86F14B857F2EB047CDF9323148BF58@prdassexch01.ad.linz.govt.nz">
      <pre class="moz-quote-pre" wrap="">Hi Even

Thanks for detailed reading and questions.

</pre>
      <blockquote type="cite">
        <pre class="moz-quote-pre" wrap="">- 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
(<a class="moz-txt-link-freetext" href="https://proj.org/operations/transformations/deformation.html">https://proj.org/operations/transformations/deformation.html</a>) 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 <a class="moz-txt-link-freetext" href="https://www.iogp.org/wp">https://www.iogp.org/wp</a>-
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.
</pre>
      </blockquote>
      <pre class="moz-quote-pre" wrap="">
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.

</pre>
      <blockquote type="cite">
        <pre class="moz-quote-pre" wrap="">- 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.
</pre>
      </blockquote>
      <pre class="moz-quote-pre" wrap="">
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 :-)

</pre>
      <blockquote type="cite">
        <pre class="moz-quote-pre" wrap="">- 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.
</pre>
      </blockquote>
      <pre class="moz-quote-pre" wrap="">
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.

</pre>
      <blockquote type="cite">
        <pre class="moz-quote-pre" wrap="">- 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.
</pre>
      </blockquote>
      <pre class="moz-quote-pre" wrap="">
Good idea - I'll build a CSV file matching the NZGD2000 deformation model.

</pre>
      <blockquote type="cite">
        <pre class="moz-quote-pre" wrap="">- 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)
</pre>
      </blockquote>
      <pre class="moz-quote-pre" wrap="">
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).

</pre>
      <blockquote type="cite">
        <pre class="moz-quote-pre" wrap="">- Any indication of the number of grids a deformation model might have ? For
the NZ case for example
</pre>
      </blockquote>
      <pre class="moz-quote-pre" wrap="">
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.

</pre>
      <blockquote type="cite">
        <pre class="moz-quote-pre" wrap="">- 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
</pre>
      </blockquote>
      <pre class="moz-quote-pre" wrap="">
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.

</pre>
      <blockquote type="cite">
        <pre class="moz-quote-pre" wrap="">- 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
</pre>
      </blockquote>
      <pre class="moz-quote-pre" wrap="">
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@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.

</pre>
      <blockquote type="cite">
        <pre class="moz-quote-pre" wrap="">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.
</pre>
      </blockquote>
      <pre class="moz-quote-pre" wrap="">
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.   <a class="moz-txt-link-freetext" href="https://www.youtube.com/watch?v=xdvI9Yrs4iA&t=141s">https://www.youtube.com/watch?v=xdvI9Yrs4iA&t=141s</a>
There is also some information on the LINZ website at <a class="moz-txt-link-freetext" href="https://www.linz.govt.nz/data/geodetic-system/datums-projections-and-heights/geodetic-datums/new-zealand-geodetic-datum-2000-nzgd2000/nzgd2000-deformation-model">https://www.linz.govt.nz/data/geodetic-system/datums-projections-and-heights/geodetic-datums/new-zealand-geodetic-datum-2000-nzgd2000/nzgd2000-deformation-model</a>,
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.

</pre>
      <blockquote type="cite">
        <pre class="moz-quote-pre" wrap="">Even
</pre>
      </blockquote>
      <pre class="moz-quote-pre" wrap="">
________________________________

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 <a class="moz-txt-link-abbreviated" href="mailto:info@linz.govt.nz">info@linz.govt.nz</a>) 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
<a class="moz-txt-link-abbreviated" href="mailto:PROJ@lists.osgeo.org">PROJ@lists.osgeo.org</a>
<a class="moz-txt-link-freetext" href="https://lists.osgeo.org/mailman/listinfo/proj">https://lists.osgeo.org/mailman/listinfo/proj</a>
</pre>
    </blockquote>
    <br>
    <pre class="moz-signature" cols="72">-- 
*************************************
Brian Shaw
Geodesist
NOAA/NOS/National Geodetic Survey
Phone # 240-533-9522</pre>
  </body>
</html>