Converting NAVD88 to EGM96 compared to vdatum tool

[tovogt]

I would like to see the difference between NAVD88 and EGM96 at a particular point (lon/lat: -80.90170, 32.03670) at the US east coast. The NOAA's vertical datum transformation tool can provide that information (with reference date 1997.0):

Now, in order to have a corresponding shell command, I would like to be able to use cs2cs for the same transformation. From my understanding, the correct command would look like this:

$ cs2cs +proj=longlat +datum=NAD83 +units=m +geoidgrids=g2018u0.gtx +to +proj=longlat +datum=WGS84 +units=m +geoidgrids=egm96_15.gtx
-80.90170 32.03670 0
-80.90	32.04 1.41

However, as you can see from the output, the value returned by proj is 1.41 while the vdatum tool gives me -0.033.

What am I doing wrong here?

[rouault]

What am I doing wrong here?

You're using the deprecated PROJ string syntax and approximate NAD83 ~= NAD83(2011) and WGS 84(G1674) ~= WGS84, which doesn't lead to the most accurate results

If you run projinfo -s "NAD83(2011) + NAVD88 height" -t "WGS84 (G1762) + EGM96 height" --spatial-test intersects, you'll get the following conversion pipeline as the priority result:

unknown id, Inverse of NAD83(2011) to NAVD88 height (3) + Conversion from NAD83(2011) (geog3D) to NAD83(2011) (geocentric) + Inverse of ITRF2008 to NAD83(2011) (1) + Inverse of WGS 84 (G1762) to ITRF2008
 (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog3D) + Inverse of WGS 84 to WGS 84 (G1762) + WGS 84 to EGM96 height (1) + WGS 84 to WGS 84 (G1762), 5.025 m, United States (USA) -
 CONUS onshore - Alabama; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Mi
chigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; Sou
th Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming., at least one grid missing

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +proj=vgridshift +grids=us_noaa_g2018u0.tif +multiplier=1
  +step +proj=cart +ellps=GRS80
  +step +inv +proj=helmert +x=0.99343 +y=-1.90331 +z=-0.52655 +rx=0.02591467
        +ry=0.00942644999999999 +rz=0.01159935 +s=0.00171504 +dx=0.00079
        +dy=-0.0006 +dz=-0.00134 +drx=6.667e-05 +dry=-0.00075744 +drz=-5.133e-05
        +ds=-0.00010201 +t_epoch=1997 +convention=coordinate_frame
  +step +inv +proj=cart +ellps=WGS84
  +step +inv +proj=vgridshift +grids=us_nga_egm96_15.tif +multiplier=1
  +step +proj=unitconvert +xy_in=rad +xy_out=deg
  +step +proj=axisswap +order=2,1

If you run
$ echo 32.03670 -80.90170 0 | PROJ_NETWORK=ON cs2cs -d 8 "NAD83(2011) + NAVD88 height" "WGS84 (G1762) + EGM96 height"
you'll get:
32.03670622 -80.90170308 -0.05035059
which is much closer to the VDatum result

As the pipeline involves a time dependent transformation between NAD83(2011) and ITRF2008, the coordinate epoch can matter (if you don't specify it, PROJ will ignore the time dependent terms, or equivalently will return as if it was 1997.0 which is the central epoch of the transformation)

If you try with 2010.0 which is the "standard" coordinate epoch for NAD83(2011), you'll get:

$ echo 32.03670 -80.90170 0 2010 | PROJ_NETWORK=ON cs2cs -d 8 "NAD83(2011) + NAVD88 height" "WGS84 (G1762) + EGM96 height"
32.03670660 -80.90170506 -0.04050265 2010

and if I try with ~ today epoch:
$ echo 32.03670 -80.90170 0 2021.8 | PROJ_NETWORK=ON src/cs2cs -d 8 "NAD83(2011) + NAVD88 height" "WGS84 (G1762) + EGM96 height"
32.03670693 -80.90170686 -0.03156374 2021.8

But I don't know if NOAA is using time dependent tranformations. You may ask them.

[tovogt]

Thanks a lot, this is really helpful!

• NOAA's vertical datum tool allows to specify the reference date. The screenshot above is with using the (default) reference date 1997.0. Is there an explanation for the remaining difference between the outputs (-0.050 from cs2cs vs. -0.033 from NOAA's tool)?
• Is there a particular reason why you use WGS84 (G1762) instead of WGS84 (G1674) which is the one used by NOAA's tool?
• Why doesn't this work equally well: cs2cs EPSG:6317+5703 EPSG:7662+5773

[rouault]

* Is there a particular reason why you use `WGS84 (G1762)` instead of `WGS84 (G1674)` which is the one used by NOAA's tool?

Ah I'm too used to use G1762. I get the same results using G1674

Why doesn't this work equally well: cs2cs EPSG:6317+5703 EPSG:7662+5773

EPSG:6317 is a geocentric CRS. You should use EPSG:6318 for the geographic 2D CRS
Same for EPSG:7662, you should use EPSG:9056

[tovogt]

Great, that works:

$ echo 32.03670 -80.90170 0 1997.0 | PROJ_NETWORK=ON cs2cs -d 8 "EPSG:6318+5703" "EPSG:9056+5773"
32.03670622	-80.90170308 -0.05035059 1997.0

Is there an explanation for the remaining difference between the outputs (-0.050 from cs2cs vs. -0.033 from NOAA's tool)? As far as I know, PROJ uses the same data as NOAA's tool for the GEOID18 model. It just converts the data to the gtx format, right?

[rouault]

Is there an explanation for the remaining difference between the outputs (-0.050 from cs2cs vs. -0.033 from NOAA's tool)?

No idea what NOAA's tool does. I've submitted the question on the mailing list on https://lists.osgeo.org/pipermail/proj/2021-October/010373.html. Got a private response from someone has NOAA who would push this to the VDatum team

PROJ >= 7 no longer uses .gtx files, but .tif files. Anyway, whatever the format, all of them should be equivalent to the original NOAA publised dataset

[gdt]

@tovogt You might consider breaking up this transform into its constituent pieces and validate each of them separately from NGS tools to proj. That might narrow down what's going on. I went down this rabbit hole trying to establish a GPS-derived NAVD88 height for a MassDOT benchmark, to compare to the published NAVD88 height determined by leveling in the 1960s, and was wicked excited to be 36mm off. I didn't deal with WGS84 geoid models though.

Starting with NAVD88, one uses GEOID18 to transform to NAD83(2011) epoch 2010.0 ellipsoidal height. I am guessing that this doesn't involve an epoch really because NAVD88 is a static CRS and NAD83(2011) epoch 2010.0 is too. I suspect NGS has an online calculator for that.

Then, there's transformation from NAD83(2011) to WGS84(G1674/G1762). My understanding is that G1674 and G1762 are within a few mm of each other and also to ITF2008, and there is no way to get accurate coordinates in WGS84 anyway. The NOAA page says that really they are using ITF2008, which makes sense, and I bet they'd say the same thing for G1762. Epoch definitely matters here and I suggest you look up nearby CORS and look at the vertical velocities in NAD83(2011) and ITRF2014 both (or IGS14, pretty much equivalently, whatever they provide for an ITRF-family frame). In New England, NAD83 velocities are about a tenth of ITRF velocities, speaking very roughly.

You might look at proj's transforms for G1762, G1674, and ITRF2008. Ideally they would all more or less be the same, but I'm unclear on that, and earlier was bitten by the NAD83/WGS84 null transform.

Finally, there's EGM96, and AIUI that was defined on WGS84(G873). But I guess you can apply it to more recent WGS84. You are asking proj and NGS for the same "add EGM96 geoid height to WGS84 HAE", so that doesn't need resolving for comparison, but it's not exactly clear what the result means. I wonder if VDATUM is trying to transform to G873 before applying EGM96. Having thought a bit I think that's the right thing to do.

[tovogt]

@gdt That's a good idea! I found that conversion from NAD83(2011) + NAVD88 height to WGS84 (G1674) yields the same result for cs2cs and NOAA's tool, but the conversion from WGS84 (G1674) to WGS84 (G1674) + EGM96 height is different, see here:

echo 32.03670622 -80.90170308 -32.97402014 | PROJ_NETWORK=ON cs2cs -d 8 "WGS84 (G1674)" "WGS84 (G1674) + EGM96 height"
32.03670622	-80.90170308 -0.05035059

What does this mean? Both proj and NOAA's tool use the same 0.25 degree reference file for the EGM96 geoid, as a tif file for proj, and as a gtx file for NOAA's tool. The values in the files differ by not more than a millimeter. We can exclude this as the source of the difference, I think.

[rouault]

We can exclude this as the source of the difference, I think.

They can possibly use a different interpolation method than the bilinear one used by PROJ. Seeing the values at the 4 corners around the point of interest, they differ by more than 60 cm, so a different interpolation method can likely result in a diffrence of a few centimers:

$  gdallocationinfo -wgs84 ../../PROJ-data/us_nga/us_nga_egm96_15.tif -80.90170308  32.03670622	
Report:
  Location: (396P,232L)
  Band 1:
    Value: -32.7865600585938
$ gdallocationinfo ../../PROJ-data/us_nga/us_nga_egm96_15.tif 397 232
Report:
  Location: (397P,232L)
  Band 1:
    Value: -33.1294937133789
$ gdallocationinfo ../../PROJ-data/us_nga/us_nga_egm96_15.tif 396 233
Report:
  Location: (396P,233L)
  Band 1:
    Value: -32.5720863342285
$ gdallocationinfo ../../PROJ-data/us_nga/us_nga_egm96_15.tif 397 233
Report:
  Location: (397P,233L)
  Band 1:
    Value: -32.9464874267578

I'd bet if you used EGM2008 instead, you might get closer results between PROJ and NOAA's tool.

[tovogt]

Okay, it's a bit unsatisfying that we don't know for sure which interpolation method is actually used by NOAA's tool, but from my side, this issue is solved. All of your comments were extremely helpful and I learned a lot about georeferencing, thanks!

[gdt]

It's an interesting question what the WGS84 defining document says about EGM96 and interpolation. Perhaps only theoretically interesting.....