[pdal] Regarding Reprojection Between Vertical Datums

[Greg Troxel]

Nicholas Stanley <nicholas.stanley at luminartech.com> writes:

> Our "ground truth" is using WGS84 Ellipsoidal Height. The survey was
> created in ~2017, but unfortunately I don't have much additional info
> about the Epoch.

One would assume that the coordinate epoch is the date of the
observations.

I am very curious about this survey and how it was done.  I do not know
of any mechanisms to obtain precise coordinates in WGS84 (assuming you
aren't working within NGA).   If it's uncorrected data I don't see how
it would be that useful.   If it's RTK or static carrier phase, I don't
see how it could be in WGS84.

> Technically, yes, I am in California - is it fair to say this
> complicates things somewhat? The data for conversion is a USGS survey
> of the San Francisco Bay Area.

As you know California is unstable :-) Most stations in the US are
located on the North American plate and their NAD83 velocities are low,
on the order of a few mm/year.  Parts of CA are on the Pacific plate and
have much higher velocities.  And, NAD83 has a velocity relative to
WGS84/ITRF.  So if you are trying to get things right at the 10 cm
level, you'll have to be careful about this.  NGS has lots of
information about seismic/plate-motion coordinate changes.  I don't
really understand the details as I'm in Massachusetts, far from a plate
boundary.  But, if you are going from NAD83(2011) epoch 2010.0 to ITRF,
I think the Pacific plate issues won't be that much trouble.

> Following the logic of your process:
>
> *Step 1: Convert the **/NAD83(2011) / Conus Albers + NAVD88 height
> /**to NAD83 Ellipsoidal using GEOID18 *
>

> E.G.: pdal translate
> USGS_LPC_CA_NoCAL_Wildfires_B5b_2018_w2274n1959.laz Step1_NAD83.laz
> reprojection --filters.reprojection.out_srs="+init=EPSG:xxxx
> +geoidgrids=~/Downloads/vdatum/core/geoid18/v2018prvi.gtx"
> --writers.las.a_srs=EPSG:xxxx

v2018prvi looks like a grid for the Puerto Rico region.   Is there a
v2018conus.gtx?  Why are you using the prvi one?

If you run pdalinfo on the .laz file, does it output a CRS, and does
that CRS match the documented CRS for the LIDAR data?

What I am also suggesting is that you obtain a point manually (some
lat/lon and the height) from your dataset, at one of the ground truth
dataset points, and then take that and do the various conversions
manually using NGS and other tools, such as
  https://www.ngs.noaa.gov/NCAT/
  https://www.ngs.noaa.gov/GEOID/GEOID18/
  https://www.ngs.noaa.gov/GEOID/GEOID18/computation.html
  https://geographiclib.sourceforge.io/cgi-bin/GeoidEval

> (Not sure which EPSG codes to use here - any guess?)

I don't think you can assume that there is an EPSG code for a compound
CRS that's e.g. a combination of a UTM zone and a height.

> *Step 2. Convert resulting Step 1. file into destination projection
> (WGS84 UTM 10 Ellipsoidal https://epsg.io/32610, or alternatively,
> https://epsg.io/7789 ITRF 2014)
>
> E.G.: pdal translate
> USGS_LPC_CA_NoCAL_Wildfires_B5b_2018_w2274n1959.laz Step1_NAD83.laz
> reprojection --filters.reprojection.out_srs="+init=EPSG:xxxx
> +geoidgrids=~/Downloads/vdatum/core/geoid18/v2018prvi.gtx"
> --writers.las.a_srs=EPSG:xxxx

I think you need to define a compound CRS that includes height and avoid
handling 3D data with a 2D CRS.

I am not really familiar with pdal translate options; your arguments
look like old-style proj rather than modern proj which lets you query
for transformations and see them.
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 194 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/pdal/attachments/20220712/0343925a/attachment.sig>