Add support for NADCON5 grids and transformation method

[rouault]

to be used with NADCON5 grids added in PROJ-data per OSGeo/PROJ-data#89. Fixes #2194

Main enhancements:

• addition of a +proj=gridshift method that can perform adjustments by applying longitude and latitude offsets (hgridshift) and/or vertical offsets (vgridshift). The and is used for NADCON5 grids which starting with NAD83(HARN) to NAD83(FBN) include ellipsoidal height offsets. The gridshift method also supports a +interpolation=bilinear/biquadratic method. By default, bilinear is applied, unless the grid has a "interpolation" metadata item set to biquadratic, which is the case for NADCON5 grids as NGS uses that method in their tools.
• modification in the EPSG importer to import transformations for the method "NADCON5 (3D)"
• refresh of the database with that
• mapping of NADCON5(2D) and NADCON5 (3D) methods to +proj=gridshift
• addition of concatenated operations do be able to do all or part of the NAD27 -> NAD83(86) -> NAD83(HARN) -> NAD83(FBN) -> NAD83(NSRS2007) -> NAD83(2011) transformation chain.

Behavior changes:

• NAD27 to NAD83 will now use the NADCON5 grids for CONUS and Alaska
• NAD83 to NAD83(HARN) will now use the NADCON5 grids instead of the HPGN grids
• NAD83 to NAD83(2011) is no longer a null transformation, but actually chains 4 successive NADCON5 grid based transformations

With all that we perfectly reproduce the results of NGS tools such as their REST API

Demo:

1. CONUS:

• NAD27 to NAD83:

$ echo 40 -80 | PROJ_DATA=/home/even/proj/PROJ-data/us_noaa:data bin/cs2cs -d 10 "NAD27" "NAD83" 
40.0000719424 -79.9997764202 0.0000000000

$ curl -s "https://geodesy.noaa.gov/api/ncat/llh?lat=40.0&lon=-80.0&inDatum=nad27&outDatum=nad83(1986)" | jq "[.destLat, .destLon, .destEht]"
[
  "40.0000719424",
  "-79.9997764202",
  "N/A"
]

• NAD83 to NAD83(2011) :

$ echo 40 -80 | PROJ_DATA=/home/even/proj/PROJ-data/us_noaa:data bin/cs2cs -d 10 "NAD83(1986)" "NAD83(2011)"
39.9999983008 -79.9999976143 0.0000000000

$ curl -s "https://geodesy.noaa.gov/api/ncat/llh?lat=40.0&lon=-80.0&inDatum=nad83(1986)&outDatum=nad83(2011)" | jq "[.destLat, .destLon, .destEht]"
[
  "39.9999983008",
  "-79.9999976143",
  "N/A"
]

• NAD83(HARN) to NAD83(2011) with ellipsoidal height offset :

$ echo 40 -80 100 | PROJ_DATA=/home/even/proj/PROJ-data/us_noaa:data bin/cs2cs -d 10 "NAD83(HARN)" "NAD83(2011)" --3d
39.9999999323 -79.9999997660 99.9755066483

$ curl -s "https://geodesy.noaa.gov/api/ncat/llh?lat=40.0&lon=-80.0&eht=100&inDatum=nad83(harn)&outDatum=nad83(2011)" | jq "[.destLat, .destLon, .destEht]"
[
  "39.9999999323",
  "-79.9999997660",
  "99.976"
]

2. Alaska:

• NAD27 to NAD83:

$ echo 61 -158 | PROJ_DATA=/home/even/proj/PROJ-data/us_noaa:data bin/cs2cs -d 10 "NAD27" "NAD83" 
60.9992880511 -158.0023022840 0.0000000000


$ curl -s "https://geodesy.noaa.gov/api/ncat/llh?lat=61.0&lon=-158.0&inDatum=nad27&outDatum=nad83(1986)" | jq "[.destLat, .destLon, .destEht]"
[
  "60.9992880511",
  "-158.0023022840",
  "N/A"
]

• NAD83(HARN) to NAD83(2011) with ellipsoidal height offset, and a grid where the ellipsoidal height offset resolution is different from the longitude&latitude offsets:

$ echo 61 -158 100 | PROJ_DATA=/home/even/proj/PROJ-data/us_noaa:data bin/cs2cs -d 10 "NAD83(1992)" "NAD83(2011)" --3d
61.0000001274 -157.9999998340 100.0632843300

$ curl -s "https://geodesy.noaa.gov/api/ncat/llh?lat=61&lon=-158&eht=100.0&inDatum=nad83(1992)&outDatum=nad83(2011)" | jq "[.destLat, .destLon, .destEht]"
[
  "61.0000001274",
  "-157.9999998340",
  "100.063"
]

[kbevers]

I think this is an amazing addition! The generalisation of gridshifting is nice. I've got a single in-lined comment. I'm not suggestion to change anything, I just need my curiosity satisfied :-)

[kbevers]

What use-case is behind this addition? A similar thing could be done with pop and push

[rouault]

The use-case is for example NAD83 to NAD83(2011) for which elliposoidal height changes are not significant given that NAD83 is 2D-only officially and the grid for the transformation NAD83 to NAD83(HARN) has no ellipsoidal height offset. Hence it is useless to pay the price for the ellipsoidal height offset interpolation for the 3 next grids of the concatenated transformation.

Below performance with the pipeline output by "projinfo -s NAD83 -t NAD83(2011)" that will insert +no_z_transform where it makes sense

$ PROJ_DATA=/home/even/proj/PROJ-data/us_noaa:data bin/bench_proj_trans \
    --pipeline "+proj=pipeline +step +proj=axisswap +order=2,1 \
                     +step +proj=unitconvert +xy_in=deg +xy_out=rad \
                     +step +proj=gridshift +grids=us_noaa_nadcon5_nad83_1986_nad83_harn_conus.tif \
                     +step +proj=gridshift +no_z_transform +grids=us_noaa_nadcon5_nad83_harn_nad83_fbn_conus.tif \
                     +step +proj=gridshift +no_z_transform +grids=us_noaa_nadcon5_nad83_fbn_nad83_2007_conus.tif \
                     +step +proj=gridshift +no_z_transform +grids=us_noaa_nadcon5_nad83_2007_nad83_2011_conus.tif \
                     +step +proj=unitconvert +xy_in=rad +xy_out=deg \
                     +step +proj=axisswap +order=2,1" 40 -100 --noise-x 1 --noise-y 1
40 -100 -> 39.9999988651031 -99.9999991253051
Duration: 3293 ms
Throughput: 1.52 million coordinates/s

compared to using push/pop and letting z interpolation:

$ PROJ_DATA=/home/even/proj/PROJ-data/us_noaa:data bin/bench_proj_trans \
       --pipeline "+proj=pipeline +step +proj=axisswap +order=2,1 \
                        +step +proj=unitconvert +xy_in=deg +xy_out=rad \
                        +step +proj=push +v_3 \
                        +step +proj=gridshift +grids=us_noaa_nadcon5_nad83_1986_nad83_harn_conus.tif \
                        +step +proj=gridshift +grids=us_noaa_nadcon5_nad83_harn_nad83_fbn_conus.tif \
                        +step +proj=gridshift +grids=us_noaa_nadcon5_nad83_fbn_nad83_2007_conus.tif \
                        +step +proj=gridshift +grids=us_noaa_nadcon5_nad83_2007_nad83_2011_conus.tif \
                        +step +proj=pop +v_3 \
                        +step +proj=unitconvert +xy_in=rad +xy_out=deg \
                        +step +proj=axisswap +order=2,1" 40 -100 --noise-x 1 --noise-y 1
40 -100 -> 39.9999988651031 -99.9999991253051
Duration: 3604 ms
Throughput: 1.39 million coordinates/s

Juts seeing that the difference is no longer that big since I added various optimizations in grid caching & interpolation, but it was more significant at an earlier stage of developpement. For some grids like Alaska where there's a dedicated subgrid for the ellipsoidal height offset, skipping the Z interpolation saves I/O, which can speed up things in a remote grid use case.

[gdt]

The use-case is for example NAD83 to NAD83(2011) for which elliposoidal height changes are not significant given that NAD83 is 2D-only officially and the grid for the transformation NAD83 to NAD83(HARN) has no ellipsoidal height offset. Hence it is useless to pay the price for the ellipsoidal height offset interpolation for the 3 next grids of the concatenated transformation.

I would expect that NAD83(1986) does not have ellipsoidal or XYZ forms. I guess you are saying that once you have started with a 2D datum the point is to skip additional work and that the user has to treat the result as 2D only even though it is a value in a 3D CRS. Did I follow that right?

[rouault]

Did I follow that right?

yes, taking into account the ellipsoidal height offsets of the HARN -> FBN -> NSRS2007 -> 2011 transformations when the starting point is 86 would give a false sense of accuracy given that the offset between 86 and HARN is probably much more bigger but nobody knows its value.

[kbevers]

Thanks for the explanation. It makes a lot of sense to optimize performance!

[gdt]

Great to see this. Not necessarily related, but I find "NAD83" confusing as sometimes it is an ensemble name and sometimes it is NAD83(1986), the original realization. That along with NAD83(1986) == WGS(TRANSIT) and similar ensemble naming issues with WGS84 leads to wrong results in older proj versions, so any steps to reduce this confusion, even if just in docs, is helpful. I lean to saying NAD83(1986) whenever the original is meant, even if that isn't what NGS does. The alternative is to say NAD83(ensemble) when that's meant. I think there is a lot of data casually labeled NAD83 that is in a later realization, and that supports the bare word being the ensemble.

Also, I think there are grids in the base proj distribution for nad27 and I wonder if those are going to be removed, if proj will no longer be using them, or if those transforms are still there but will be lower priority unless explicitly used.

[rouault]

I think there is a lot of data casually labeled NAD83 that is in a later realization, and that supports the bare word being the ensemble.

I'm not in position to solve this. The best we can do is adhere with what official authorities say. Otherwise we'd add even more confusion to something confusing.

Also, I think there are grids in the base proj distribution for nad27 and I wonder if those are going to be removed, if proj will no longer be using them, or if those transforms are still there but will be lower priority unless explicitly used.

The old grids and corresponding database records will still be there. They would be used if the NADCON5 grids are not available locally and PROJ networking is not enabled, or if the user specifically selects those transformations in the list of transformations returned by our API that return a list of transformations, or if the user specifies an explicit pipeline that use them.

[gdt]

I know you can't solve it either; I was more suggesting that proj within the project be consistent and adopt two names, one for the ensemble and one for NAD83(1986).

[busstoptaktik]

I was more suggesting that proj within the project be consistent and adopt two names, one for the ensemble and one for NAD83(1986)

@gdt that would probably be hard to do in a consistent way, since the de facto convention of postfixing a paranthesized realization epoch postdates NAD83, so NAD83(1986) is probably an afterthought, first being meaningful after the introduction of NAD83(HARN). Perhaps, the NSRS really should be considered the ensemble name, and the various NAD83 realizations considered NSRS-realizations? (yes - that was a question: I'm not very well versed in US reference frames and systems)

[rouault]

NAD83, so NAD83(1986)

Note: in the EPSG dataset, "NAD83(1986)" is an official alias to "NAD83" (the later is the official/main name) for both the CRS and the datum objects. You can use both names indifferently when using PROJ functions that accept names as inputs.

Regarding ensembles, I doubt this is going to change, as it would be only for those past CRS, whereas the focus of NGS is now oriented towards the future NATRF2022 and associated datums, planned for 2025, where there will be I believe a ~ 1.5m offset in coordinates (due to the "incorrect" centre of the Earth used for the NAD83 family), so it is unlikely that the 2022 systems could be part of an hypothetical NAD83 ensemble.

[rouault]

Merging