Vertical datums removed in .ept files?

[mattbeckley]

I was curious as to how vertical datums were handled in the ept data. I understand that all the ept has been re-projected into web mercator. But, some of the original data that I have been looking at in:
ftp://rockyftp.cr.usgs.gov/vdelivery/Datasets/Staged/Elevation/LPC/Projects/
have geoids applied (GEOID09, GEOID06, etc). Also, not all of them use meters for the units.

My question is: When the data were translated to EPSG:3857 for the ingest into .ept tiles, was the vertical datum accounted for (i.e. was the height converted back to an ellipsoidal height as opposed to an orthometric height?)

Any info or insight would be appreciated. Thanks!

[hobu]

When the data were translated to EPSG:3857 for the ingest into .ept tiles, was the vertical datum accounted for (i.e. was the height converted back to an ellipsoidal height as opposed to an orthometric height?)

Sorry for the trouble on this. Our intention was to normalize all of the vertical to ellipsoid as part of our processing. We had an issue with PDAL/GDAL/PROJ doing that and we backed off the normalization of the verticals during the final processing of the collects to EPT. To answer your question, the heights should be whatever orthometric heights they were provided in, and the metadata for the VERTCS should be available in the ept-sources of the EPT object.

We will try to script up a process to pull up those VERTCS entries into the top-level EPT metadata. If there are mismatches in any of the sources, we're in a bit of a bind though 😦 Most of the collections that we didn't process to EPT were left alone for SRS metadata consistency issues, so hopefully the challenging ones are left to that current set, not the set of collects that have been processed to EPT.

[mattbeckley]

Hi Howard. Thanks for the response, and I understand that it would be difficult to normalize all the vertical datums, but I just wanted to double check. I pulled all the json files with the ept-sources directory for each dataset, and checked the vertical datums. There was only 1 dataset that had some mixed datums. OK_GrandLake_2010 has some datasets using GEOID03, and some using GEOID09
['VERT_CS["NAVD88 - Geoid09 (Feet)"',
['VERT_CS["NAVD88 - Geoid03 (Feet)"']

But thankfully, that was the only one that my script found. Thanks for your work on this project, it's a great dataset!

[hobu]

Can you post the script or a JSON map of collect=>vertcs? We'll use that information to go through and update the resources.

I think we should remove OK_GrandLake_2010.

[mattbeckley]

Attached is the python script I used to loop through and pull out the vertical CRS info. some of the JSON files had some non-ascii characters that were messing things up for me, so also attached is the bash script that fixes that. I'm not a software engineer by any means, so keep that in mind :)

GetVCRSInfo.zip

[mattbeckley]

In revisiting this issue, I had another question. I understand that the elevations retained their original vertical datums, but were the units for the elevations all converted to meters? I am looking at a dataset:
USGS_LPC_VA_ChesapeakeBayNorth_2015_LAS_2017
the original data from USGS is in NAD83(2011) / Virginia North (ftUS) with NAVD88 vertical datum applied, with both in units of US survey feet. This is also what the metadata shows in the ept-sources files in the AWS bucket. However, when I download data from the AWS bucket:

{
  "pipeline": [
    {
      "type": "readers.ept",
        "filename": "https://s3-us-west-2.amazonaws.com/usgs-lidar public/USGS_LPC_VA_ChesapeakeBayNorth_2015_LAS_2017",
      "bounds": "([-8769205, -8765318], [4634035,4636958])"
    },
     "test_VA.laz"
  ]
}

, the Z units appear to be in meters. So, I was just curious if all elevations were converted to meters as part of the conversion to ept format?

[hobu]

the Z units appear to be in meters. So, I was just curious if all elevations were converted to meters as part of the conversion to ept format?

It was unintended, but the answer to this is yes. GDAL's units for EPSG:3857 is meters and it assumed that for vertical too and applied the math. I'm not sure I like that ☹️

[dshean]

Hey @hobu. We are revisiting this issue, and attempting to define a PROJ pipeline and 3D transformation to go from the 3DEP EPT json to local UTM zone relative to WGS84 (G2139) (https://github.com/uw-cryo/3D_CRS_Transformation_Resources/blob/main/examples/UTM_10N_WGS84_G2139_3D.wkt).

We get the following when we run:
projinfo -s EPSG:3857+5703 -t "$(cat UTM_10N_WGS84_G2139_3D.wkt)" -o PROJ --hide-ballpark --spatial-test intersects

Operation No. 1:

unknown id, Inverse of Popular Visualisation Pseudo-Mercator + Inverse of NAD83(2011) to WGS 84 (1) + Inverse of NAD83(2011) to NAVD88 height (3) + Conversion from NAD83(2011) (geog3D) to NAD83(2011) (geocen
tric) + Inverse of ITRF2014 to NAD83(2011) (1) + Inverse of WGS 84 (G2139) to ITRF2014 (1) + Conversion from WGS 84 (G2139) (geocentric) to WGS 84 (G2139) (geog3D) + UTM zone 10N, 2.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; Mic
higan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Ca
rolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming.

PROJ string:
+proj=pipeline
+step +inv +proj=webmerc +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84
+step +proj=vgridshift +grids=us_noaa_g2018u0.tif +multiplier=1
+step +proj=cart +ellps=GRS80
+step +inv +proj=helmert +x=1.0053 +y=-1.90921 +z=-0.54157 +rx=0.02678138
+ry=-0.00042027 +rz=0.01093206 +s=0.00036891 +dx=0.00079 +dy=-0.0006
+dz=-0.00144 +drx=6.667e-05 +dry=-0.00075744 +drz=-5.133e-05
+ds=-7.201e-05 +t_epoch=2010 +convention=coordinate_frame
+step +inv +proj=cart +ellps=WGS84
+step +proj=utm +zone=10 +ellps=WGS84

I'm surprised by the 2.025 meter uncertainty (presumably from the EPSG:3857 definition), I'm not sure we need the helmert here, as I think this was done when you went from the original laz (NAD83(2011)+NAVD88) to the web mercator horizontal CRS? But the vertical CRS of the EPT should still be NAVD88 (presumably geoid2018)?

Some of the other candidate PROJ operations did not include the helmert, but they reverted to using older NAVD88 grids, e.g.:

Operation No. 2:

unknown id, Inverse of Popular Visualisation Pseudo-Mercator + Inverse of NAD83(NSRS2007) to WGS 84 (1) + Inverse of NAD83(NSRS2007) to NAVD88 height (1) + NAD83(NSRS2007) to WGS 84 (1) + WGS 84 to WGS 84 (G2139) + UTM zone 10N, 6.05 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; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming.

PROJ string:
+proj=pipeline
+step +inv +proj=webmerc +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84
+step +proj=vgridshift +grids=us_noaa_geoid09_conus.tif +multiplier=1
+step +proj=utm +zone=10 +ellps=WGS84

When we test the conversion...
echo 8274397.455 -19995708.608 1400 | PROJ_NETWORK=ON cs2cs -f "%.3f" -r EPSG:3857+5703 "$(cat UTM_10N_WGS84_G2139_3D.wkt)"

Using coordinate operation Inverse of Popular Visualisation Pseudo-Mercator + WGS 84 to WGS 84 (G2139) + Transformation from NAVD88 height to WGS 84 (ballpark vertical transformation, without ellipsoid height to vertical height correction) + UTM zone 10N
-2397786.144 7990294.726 1400.104

So not actually using that first proj pipeline operation, and not applying the ~20 m vertical shift required for this region.

Hoping you or others can provide some guidance here. We are considering reverting to the OPR laz or the LPC laz available on the old ftp server, so we can get this right, and avoid uncertainty from intermediate transformations. Thanks!
@kvenkman

[dshean]

We are using the WA_MtBaker_2015 dataset for these tests:
https://s3-us-west-2.amazonaws.com/usgs-lidar-public/USGS_LPC_WA_MtBaker_2015_LAS_2017/ept.json
https://rockyweb.usgs.gov/vdelivery/Datasets/Staged/Elevation/LPC/Projects/WA_MtBaker_2015/WA_MtBaker_2015/LAZ/USGS_LPC_WA_MtBaker_2015_10UEU8198.laz

[dshean]

Hey @hobu (or @mattbeckley, @cjcrosby), any insight? We're kind of in a holding pattern right now, and at a design turning point - need to decide whether we can proceed with the ept files, or we need to go back to original laz tiles to avoid CRS ambiguity issues. Thanks!

[hobu]

go back to original laz tiles to avoid CRS ambiguity issues

The CRS ambiguity started at the original LAZ tiles 😄

[hobu]

The verticals on the EPT data should be untouched during the conversion process to 3857. That said, they are inconsistently defined on those original LAZ tiles too. At one time, I made Plumbline to go back to the SRTM data and do some simple comparison of points to try to guess what the vertical of various data might be.

It would be handy to have a key to map all of the resources and the vertical CRS information so we could put it back into the EPT data. Money/time/effort thing, of course.

[dshean]

OK, thanks @hobu. I agree that would be useful, as would some additional sanity checks and intelligent guesses (maybe against COP30 this time around?).

I guess I'm still not clear on whether your transformation to EPSG:3857 included the NAD83(2011) to WGS84 ellipsoid helmert transformation (+x=1.0053 +y=-1.90921 +z=-0.54157...), which involves an additional vertical height offset. Do we need to include this helmert transformation, or do we just need to simply remove the NAVD88 offset from the EPT points? I'm guessing that there was no helmert if you used EPSG:3857 WGS84 Ensemble and not a specific realization. Maybe you can point me to the code you used if public, or check the logs?

I'm also wary of gotchas involving the spherical/ellipsoidal mix for EPSG:3857. Maybe a non-issue as underlying horizontal CRS is WGS84 Ellipsoid.

Thanks!

[hobu]

The data are processed using the connormanning/entwine docker container, so that is the place to look for information about PROJ, PDAL, and Entwine versions to try to decipher what's going on.