[pdal] Regarding Reprojection Between Vertical Datums

Nicholas Stanley nicholas.stanley at luminartech.com
Tue Jul 12 22:32:59 PDT 2022 

*
For the fun of it, here's a link to the Ground Truth and Data to be 
reprojected in alignment to it. If Someone thinks they can crack this, 
please give it a shot!

https://drive.google.com/drive/folders/1KdUc1xnOtNl0RkWfr1-C_gh7jzU7q9U1?usp=sharing


*

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


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",*

*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":
     {
*  "EPSG:4326":*
       {
         "bbox":
         {
           "maxx": -122.4022284,
           "maxy": 37.81572701,
           "maxz": 99.83,
           "minx": -122.4164013,
           "miny": 37.80469945,
           "minz": -70.28
         },
         "boundary": { "type": "Polygon", "coordinates": [ [ [ 
-122.413254958928221, 37.804699453399095, -70.28 ], [ 
-122.416401346373235, 37.813279215545926, -70.28 ], [ 
-122.405373596515062, 37.815727007213006, 99.83 ], [ 
-122.402228384472423, 37.807146965417473, 99.83 ], [ 
-122.413254958928221, 37.804699453399095, -70.28 ] ] ] }
       },
       "native":
       {
         "bbox":
         {
           "maxx": -2273000,
           "maxy": 1959999.99,
           "maxz": 99.83,
           "minx": -2273999.99,
           "miny": 1959000,
           "minz": -70.28
         },
         "boundary": { "type": "Polygon", "coordinates": [ [ [ 
-2273999.990000000223517, 1959000.0, -70.28 ], [ 
-2273999.990000000223517, 1959999.99, -70.28 ], [ -2273000.0, 
1959999.99, 99.83 ], [ -2273000.0, 1959000.0, 99.83 ], [ 
-2273999.990000000223517, 1959000.0, -70.28 ] ] ] }
       }
     },


*EPSG 4326 is a WGS84 code, but this is merely defining the bounding box 
container it seems.

*
*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"]]]
variable length header record 4 of 4:
*
---------

///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

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



On 7/12/22 16:35, Greg Troxel wrote:
> 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 Ellipsoidalhttps://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 --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/pdal/attachments/20220712/360d1b10/attachment-0001.htm>

More information about the pdal mailing list