[pdal] Regarding Reprojection Between Vertical Datums

[Greg Troxel]

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

> *The Source Materials (to be corrected) are described here:*
> https://www.fisheries.noaa.gov/inport/item/59131

Two interesting things:

  - It says vertical is NAVD88, but it also says GEOID12B.  That's a
    clue that they did their control in NAD83 lat/lon/HAE and then
    transformed to NAVD88 via GEOID12B (as one would expect given the
    publication date of the dataset).  So if you are trying to transform
    back to HAE, you should use GEOID12B rather than GEOID18.  But I
    suspect the differences between 12/18 are mostly less than 10cm, and
    they certainly won't be 2m.

  - They did ground control to reference the data, and then check
    points.  Their data quality statement claims RMS vertical errors
    better than 10 cm.

As a data point, the Massachusetts state GIS (MassGIS) has orthoimagery
(captured in conjunction with USGS), and has a similar control statement
about check points and RMS horizontal errors.  My own observations are
consistent with the published accuracy statement.  So I think it's very
likely this data set meets the quality standards it claims to.

It sounds like you have some additional data points of lat/lon/HAE or
lat/lon/orthometric-height, collected somehow, that you are calling
"ground truth", and you are trying to compare to verify the LIDAR
imagery.  You haven't explained what the collection methodology was, and
you haven't answered my question about how it can be both accurate and
in "WGS84".

> Coordinate Reference System
> CRS Type:    Geographic 3D
> *EPSG Code:    EPSG:6319*
> EPSG Name:    NAD83(2011)
> Datum Name:    NAD83 (National Spatial Reference System 2011)
> Ellipsoid:    GRS 1980
> Semimajor Axis:    6378137.0
> Flattening Ratio:    298.257222101
> *"ldrgeoid" : "GEOID 12B",*

So that says 6319 which implies HAE and then gives a GEOID, which means
it isn't really 6319 but really EPSG:6349.  And it's projected, which
means the horizontal is EPSG:6350.  I am unable to find an EPSG code for
6350+NAVD88.

Somehow, you have to ensure that pdal uses the right CRS, and that may
involve defining your own compound CRS and specifying it. 

> *That said, when i run pdal info on the file, i see the following:*
>
> "file_size": 10024109,
>   "filename": "USGS_LPC_CA_NoCAL_Wildfires_B5b_2018_w2274n1959.laz",
>   "now": "2022-07-12T21:37:52-0700",
>   "pdal_version": "2.2.0 (git-version: Release)",
>   "reader": "readers.las",
>   "stats":
>   {
>     "bbox":
>     },

bbox elided, but there is no CRS for the file shown.  I don't know why,
and where the problem is.  (I acknowledge the comment about bbox being
in WGS84 and native both, and wonder if that's in the file or from
pdalinfo.)

> *Lasinfo shows:
>
>   description          'OGC Transformation Record'
>     WKT OGC COORDINATE SYSTEM:
>     COMPD_CS["NAD83(2011) / Conus Albers + NAVD88
> height",PROJCS["NAD83(2011) / Conus
> Albers",GEOGCS["NAD83(2011)",DATUM["NAD83 (National Spatial Reference
> System 2011)",SPHEROID["GRS
> 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","1116"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","6318"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",29.5],PARAMETER["standard_parallel_2",45.5],PARAMETER["latitude_of_center",23],PARAMETER["longitude_of_center",-96],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["meter",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","6350"]],VERT_CS["NAVD88
> height",VERT_DATUM["North American Vertical Datum
> 1988",2005,AUTHORITY["EPSG","5103"]],UNIT["meter",1,AUTHORITY["EPSG","9001"]],AXIS["Up",UP],AUTHORITY["EPSG","5703"]]]

Well, that is the compound CRS I was talking about!
It looks like pdalinfo isn't printing the CRS.

> ///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///
>
> //
>
> Having done this, I see coordinates that more closely align in
> Elevation to our ground Truth (although they still appear to be off by
> ~2-3 
> meters, not 30)
>
>    Output from GEOID18
>
>                                          latitude        longitude       N         error (95% confidence interval)
> Station Name                            ddd mm ss.sssss ddd mm ss.sssss  meters    meters
> USER LOCATION                            37 48 25.00920 237 35 38.63040  -32.486    0.031

That output is in NAD83(2011) epoch 2010.0 lat/lon/HAE.

> Our "Ground Truth" in this area shows an Elevation at -29.17m.

But you said that your "Ground Truth" is in WGS84.  WGS84 HAE and NAD83
HAE are different, and you didn't transform from one to the other.  3m
seems high for a delta there, but 2m wouldn't surprise me.

-------------- 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/20220713/54bb2e3f/attachment.sig>