[PROJ] geoid model + Helmert

[Greg Troxel]

"Lesparre, Jochem" <Jochem.Lesparre at kadaster.nl> writes:

> Even Rouault wrote:
>> yes, your analysis is fully correct , and your enhanced pipeline too. However that would suck a bit to have to implement it like that as
>> doing and undoing the vgridshift is inefficient. The ideal solution would use some storage space as we discussed a few months ago but
>> couldn't reach a conclusion how to implement that.  For 3D operations only like that, we could use the v_4 component as the temporary
>> storage using push +v_3 / axisswap order=3,4 / pop +v_3 trickery
>
> Greg Troxel wrote:
>> What I think you want is to get a grid file that is indexed in CH1903+ and that gives the difference from CH1903+ HAE to LHN95 height (or
>> defines LHN95, giving an estimate from CH1903+ HAE to LN02).   I don't
>> understand why that isn't what's available.
>
> Probably that isn't available because the usual way to combine a geoid
> model and a Helmert transformation is to base them both on a
> international/regional CRS (e.g. ETRS89). I think that this problem
> should be fixed in the PROJ code, not by asking national authorities
> to make redundant geoid models.

In the US (what I've mostly paid attention to) we have a regional datum
but more or less no local datums, or really that the older datums aren't
in use much (NAD27/NGVD29).  At least I am not aware of geoid models for
NGVD29.  I'm viewing NAD83(2011/CSRS-recent) and ETRS89 as occupying
similar logical positions.

> I'm working for the national authority of the Netherlands (NSGI). We
> find using the v_4 component as the temporary storage too confusing,
> when also publishing time-dependent transformations. Therefore, we
> publish a pipeline with the inefficient doing and undoing the
> vgridshift. However, I'm still hoping for a conclusion for the
> discussion of a few months ago on how to implement the ideal solution.

I am having trouble understanding, so a brief question if you're
willing.

  I think we're talking about
   ch_swisstopo_chgeo2004_ETRS89_LHN95.tif
  which is the same data aas
   https://opendata.swiss/en/dataset/geoidmodell-chgeo2004-in-etrs89
  which says
    Swiss geoid model in the version of 2004 (CHGeo2004) in the reference
    system ETRS89. It forms the zero reference surface (approximated mean
    sea level) for height determination and serves for the transformation
    between ellipsoidal heights and orthometric heights (LHN95). The geoid
    heights are stored in a regular 30 x 30 seconds grid and are
    interpolated by the biquadratic method. The model is available in
    various reference systems and data formats.

  The swisstopo grid shift file is indexed horizontally in ETRS89.
  That's easy to understand.

  The values are expressed in ETRS89 HAE, and are the ellipsoidal
  heights of the LHN95 zero surface.

So, if you want to start in ETRS89, and end up with horizontal in
CH1903+, and vertical in LHN95, you have to:

  A) Use the geoid model to get LHN95 height from the ETRS89 HAE, based on
  ETRS89 horizontal.

  B) Take the original ETRS89 LLh and Helmert transform to CH1903+, and
  then discard the h, and swap in the LHN95 H from step A.

And I think that

  C) Take ETRS89 LLh, and apply the geoid model to get LLH in
  ETRS89+LHN95, and then Helmert transform that LLH as if it were ETRS89
  LLh to get CH1903+, and treat it as CH1903+ + LHN95, then you'll end
  up with an introduced error which is the difference in ETRS89 and
  CH1903+ ellipsoid heights.

What I'm having trouble with is that for a station, I expect ETRS89 HAE
and CH1903+ HAE to be different, and that therefore one has to do the
A/B procedure above.   I don't understand how that relates to "Undoing"
the grid shift.  Perhaps these are two ways of expressing the same thing
and that's what the push/pop is all about.


I can certainly see your point that proj should work well with the grid
files as published.