[pdal] Regarding Reprojection Between Vertical Datums
Kirk Waters - NOAA Federal
kirk.waters at noaa.gov
Wed Jul 13 05:16:48 PDT 2022
Nicholas,
Based on your link to the metadata, it looks like you have at least found a
reference to the data on the NOAA system. That link describes the data as
it exists in the backend of the Data Access Viewer
<https://coast.noaa.gov/dataviewer/#/lidar/search/where:ID=9067> (the link
will zoom to that data set). But, if your original data is in Albers, it
sounds like you picked up the data from USGS. The quick cheat is to use the
Data Access Viewer and request the data in WGS84/ITRF ellipsoid heights for
the area you care about.
I think your primary issue is that the geoid model isn't getting removed.
As mentioned, you have to specify the grid to use in PDAL. The CONUS area
is split into 8 grids as distributed, though it's possible that someone has
combined them. You'll need to make sure you specify the right grid. Looking
at the metadata you linked, I see that we removed GEOID12B, so you need to
make sure that's the grid set you use (not GEOID18). The gird name will
likely look like g2012buX.gtx where X is a number from 1 to 8. The u is
what is indicating that it is CONUS. PDAL may allow you to specify multiple
grids at a time so you don't have to figure out which of the 8 it is. I
suspect it is either grid 1 or 2 for California. I noticed that VDatum also
has a g2012bu0.gtx which appears to be all 8 grids combined, so that would
be an easy pick.
Once your sure you've got the geoid taken care of, I think there are still
some other issues. As was mentioned by others, you need to specify that
you're going to transform to ITRF2008 or ITRF2014 instead of simply WGS84.
That's because PROJ is using a ballpark offset to convert to WGS84 (you can
see that with projinfo -s EPSG:6318 -t EPSG:4326) with an accuracy around 2
meters. I suspect it's using the NAD83(86) to WGS84(Transit), but doesn't
have a way to convert NAD83(2011) to NAD83(86) because the only published
way is a multi-grid hop. Something like
NAD83(2011)->NAD83(2007)->NAD83(HARN)->NAD83(86), but I might have missed a
hop in there. However, there are published Helmert transforms from
NAD83(2011) to ITRF(2008) and ITRF(2014). The accuracy of those transforms
will be much better.
One last issue. It can be hard to get the georeferencing right for
ellipsoid heights. There is no vertical code in EPSG for it. There isn't,
AFAIK, either a convenient compound CRS EPSG code or a combination of
horizontal and vertical codes to give for what you want as your output. You
can get to geographic coordinates and ellipsoid heights using the 3D codes
(e.g. for ITRF2014 it's EPSG:7912), but not projected coordinates. You may
need to specify your transform in pdal to go to 7912 and then transform to
UTM while avoiding a vertical transformation. The file georeferencing on
your output may be missing the vertical, though I haven't tried it.
Kirk Waters, PhD
NOAA Office for Coastal Management
Applied Sciences Program
Phone: 843-284-6962 (New as of 2022)
coast.noaa.gov/digitalcoast
Thunder is good, thunder is impressive; but it is lightning that does the
work - Mark Twain
On Wed, Jul 13, 2022 at 1:33 AM Nicholas Stanley <
nicholas.stanley at luminartech.com> wrote:
>
>
>
>
>
>
> * 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
> <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> <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.
>
> _______________________________________________
> pdal mailing list
> pdal at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/pdal
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/pdal/attachments/20220713/2042dc4f/attachment-0001.htm>
More information about the pdal
mailing list