Raw LAS WKT seem not to render with projinfo

[rsbivand]

In edzer/sp#75, the reported issue shifted from a problem between two R packages which is not relevant here to one of ingesting LAS 1.4 WKT strings. Comment:

None of these raw strings from LAS files (1.4, WKT specified) seem to render as expected in projinfo:

COMPD_CS["Projected", PROJCS["UTM_10N", GEOGCS [ "WGS84", DATUM [ "WGS84", SPHEROID ["WGS 84", 6378137.000, 298.257223563 ], TOWGS84 [ 0.000, 0.000, 0.000, 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000 ] ], PRIMEM [ "Greenwich", 0.000000 ], UNIT [ "metres", 1.00000000] ], PROJECTION["Transverse_Mercator"], PARAMETER["Latitude_of_Origin",0.0000000000], PARAMETER["Central_Meridian",-123.0000000000], PARAMETER["Scale_Factor",0.9996000000], PARAMETER["False_Easting",500000.000], PARAMETER["False_Northing",0.000], UNIT [ "metres", 1.00000000]] ], VERT_CS["NAVD88 (Geoid03) ContUS", VERT_DATUM["./Resources/CoordSysData/navd88_geo03_contus.bin", 1 ], UNIT [ "metres", 1.00000000] ] ] 

PROJCS["NAD83 (2011) / Conus Albers",GEOGCS["GCS_NAD_1983_2011",DATUM["D_NAD_1983_2011",SPHEROID["GRS_1980",6378137.0,298.257222101,AUTHORITY["EPSG",7019]],AUTHORITY["EPSG",1116]],PRIMEM["Greenwich",0.0,AUTHORITY["EPSG",8901]],UNIT["Degree",0.0174532925199433,AUTHORITY["EPSG",9102]],AUTHORITY["EPSG",6318]],PROJECTION["Albers",AUTHORITY["Esri",43007]],PARAMETER["False_Easting",0.0,AUTHORITY["Esri",100001]],PARAMETER["False_Northing",0.0,AUTHORITY["Esri",100002]],PARAMETER["Central_Meridian",-96.0,AUTHORITY["Esri",100010]],PARAMETER["Standard_Parallel_1",29.5,AUTHORITY["Esri",100025]],PARAMETER["Standard_Parallel_2",45.5,AUTHORITY["Esri",100026]],PARAMETER["Latitude_Of_Origin",23.0,AUTHORITY["Esri",100021]],UNIT["Meter",1.0,AUTHORITY["EPSG",9001]]],VERTCS["NAVD_1988",VDATUM["North_American_Vertical_Datum_1988",AUTHORITY["EPSG",5103]],PARAMETER["Vertical_Shift",0.0,AUTHORITY["Esri",100006]],PARAMETER["Direction",1.0,AUTHORITY["Esri",100007]],UNIT["Meter",1.0,AUTHORITY["EPSG",9001]],AUTHORITY["EPSG",5703]] 

shows two WKT strings extracted from headers of LAS files in the wild, discussed in the sp issue. I cannot find a way of running them through projinfo to try to see what gets included; the arbitrary use/omission of the vertical CRS seems worrying. Has the use of LAS WKT been examined, and what should I have been doing?

[rsbivand]

Using GDAL 3.0.4 and PROJ 7.0.0, and rgdal::showSRID(), I see the vertical CRS of the Conus Albers case discarded, and in the UTM_10N case, I see importFromWkt() failing (latest rgdal R-Forge SVN commit). The steps in R call GDAL's importFromWkt() - is it possible that proj_create() might behave differently?

Is it possible that the misbehaviour of projinfo is caused by string length?

[rsbivand]

Maybe the Conus Albers in wkt1 can be read by proj_create():

> list_coordOps(wkt1, "EPSG:4326")
Candidate coordinate operations found:  7 
Strict containment:  FALSE 
Visualization order:  TRUE 
Source: PROJCS["NAD83 (2011) / Conus Albers",
        GEOGCS["GCS_NAD_1983_2011", DATUM["D_NAD_1983_2011",
        SPHEROID["GRS_1980", 6378137.0, 298.257222101,
        AUTHORITY["EPSG", 7019]], AUTHORITY["EPSG", 1116]],
        PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG", 8901]],
        UNIT["Degree", 0.0174532925199433, AUTHORITY["EPSG",
        9102]], AUTHORITY["EPSG", 6318]], PROJECTION["Albers",
        AUTHORITY["Esri", 43007]], PARAMETER["False_Easting",
        0.0, AUTHORITY["Esri", 100001]],
        PARAMETER["False_Northing", 0.0, AUTHORITY["Esri",
        100002]], PARAMETER["Central_Meridian", -96.0,
        AUTHORITY["Esri", 100010]],
        PARAMETER["Standard_Parallel_1", 29.5,
        AUTHORITY["Esri", 100025]],
        PARAMETER["Standard_Parallel_2", 45.5,
        AUTHORITY["Esri", 100026]],
        PARAMETER["Latitude_Of_Origin", 23.0, AUTHORITY["Esri",
        100021]], UNIT["Meter", 1.0, AUTHORITY["EPSG", 9001]]],
        VERTCS["NAVD_1988",
        VDATUM["North_American_Vertical_Datum_1988",
        AUTHORITY["EPSG", 5103]], PARAMETER["Vertical_Shift",
        0.0, AUTHORITY["Esri", 100006]], PARAMETER["Direction",
        1.0, AUTHORITY["Esri", 100007]], UNIT["Meter", 1.0,
        AUTHORITY["EPSG", 9001]], AUTHORITY["EPSG", 5703]]
Target: EPSG:4326 
Best instantiable operation has accuracy: 2 m
Description: Inverse of unnamed + Inverse of NAD83 to NAD83(2011) (1) +
             NAD83 to WGS 84 (18) + axis order change (2D)
Definition:  +proj=pipeline +step +inv +proj=aea +lat_0=23 +lon_0=-96
             +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +ellps=GRS80
             +step +proj=hgridshift +grids=us_noaa_FL.tif +step
             +proj=unitconvert +xy_in=rad +xy_out=deg
Operation 3 is lacking 1 grid with accuracy 2 m
Missing grid: us_noaa_ethpgn.tif 
URL: https://cdn.proj.org/us_noaa_ethpgn.tif 
Operation 4 is lacking 1 grid with accuracy 2 m
Missing grid: us_noaa_lahpgn.tif 
URL: https://cdn.proj.org/us_noaa_lahpgn.tif 
Operation 5 is lacking 1 grid with accuracy 2 m
Missing grid: us_noaa_mshpgn.tif 
URL: https://cdn.proj.org/us_noaa_mshpgn.tif 
Operation 6 is lacking 1 grid with accuracy 2 m
Missing grid: us_noaa_alhpgn.tif 
URL: https://cdn.proj.org/us_noaa_alhpgn.tif 

Subsetting the input string to drop the COMP...[] bracketting, I see:

> list_coordOps(substring(wkt, 23, nchar(wkt)-1), "EPSG:4326")
Candidate coordinate operations found:  1 
Strict containment:  FALSE 
Visualization order:  TRUE 
Source: PROJCS["UTM_10N", GEOGCS [ "WGS84", DATUM [ "WGS84", SPHEROID
        ["WGS 84", 6378137.000, 298.257223563 ], TOWGS84 [
        0.000, 0.000, 0.000, 0.0000000000, 0.0000000000,
        0.0000000000, 0.0000000000 ] ], PRIMEM [ "Greenwich",
        0.000000 ], UNIT [ "metres", 1.00000000] ],
        PROJECTION["Transverse_Mercator"],
        PARAMETER["Latitude_of_Origin", 0.0000000000],
        PARAMETER["Central_Meridian", -123.0000000000],
        PARAMETER["Scale_Factor", 0.9996000000],
        PARAMETER["False_Easting", 500000.000],
        PARAMETER["False_Northing", 0.000], UNIT [ "metres",
        1.00000000]] ], VERT_CS["NAVD88 (Geoid03) ContUS",
        VERT_DATUM["./Resources/CoordSysData/navd88_geo03_contus.bin",
        1 ], UNIT [ "metres", 1.00000000] ]
Target: EPSG:4326 
Best instantiable operation has only ballpark accuracy 
Description: Inverse of unnamed + Transformation from WGS84 to WGS84 + axis
             order change (2D)
Definition:  +proj=pipeline +step +inv +proj=tmerc +lat_0=0
             +lon_0=-7047.38088010913 +k=0.9996 +x_0=500000
             +y_0=0 +ellps=WGS84 +step +proj=unitconvert
             +xy_in=rad +xy_out=deg

I can see something that ballpark-works, because the LAS WKT includes an unavailable grid referred to as local file in the input.

Are there guidelines for working with LAS WKT input strings? R packages using LAS are using R-spatial packages sp and sf too, and these are using PROJ and GDAL. Should I raise this on the PROJ mailing list?

[Jean-Romain]

For more information about LAS WKT see the format specification section 3.2 here.

The specification states (among other)

For definition of WKT, we refer to Open Geospatial Consortium (OGC) specification “OpenGIScoordinate transformation service implementation specification” revision 1.00 released 12 January2001, section 7 (“Coordinate Transformation Service Spec"). As there are a few dialects of WKT,please note that LAS is not using the “ESRI WKT” dialect, which does not include TOWGS84 andAuthority nodes.

But in practice it is a string so there is no technical guard and a data provider could put a string that is not compliant with the specs.

[rouault]

Those WKT strings aren't compliant with the http://portal.opengeospatial.org/files/?artifact_id=999 spec (nor any other WKT spec I'm aware of)

• The first string has the following structure COMPD_CS["bla", PROJCS[...] ], VERT_CS[...] ]. Note the imbalance of '[' and ']'. There's an extra ']', either after PROJCS[] or after VERT_CS[]. The result of this is that the VERT_CS is seen at the top-level, and not a child of the COMPD_CS. If the ']' just after PROJCS[...] was removed, then PROJ can makes sense of this WKT. I'm not super enthusiastic about supporting this badly defined WKT because it will complicate the WKT parser, unless it is proven that there's a fair number of such invalid WKT in the wild, and not just one.

• The second one has the structure PROJCS[], VERTCS[]. Which is also illegal WKT. One side issue here is the use of VERTCS instead of VERT_CS. That construct could probably be recognized more easily. But same consideration here: is it worth adding a hack for that?

[rouault]

CC @hobu who has likely seen his share of crazy LAS WKT...

[Jean-Romain]

I think this answers the question. There are also many LAS file in the wild that record non-existing epsg code so I'm not surprised.

[rsbivand]

Thanks, GDAL's SRS importFromWkt() code found the unbalanced bracket, and stripping the COMPD_CS[] also helped. Is there other WKT debugging code say in PROJ's proj_create()? The use of rgdal::list_coordOps() presents the strings to proj_create(), which seems to get further than GDAL's importFromWkt() once the strings are tidied. Changing to VERTCS from VERT_CS doesn't help since the grid file is local.

So maybe there is a way of warning users that their LAS coordinate reference system fails on coercion to sp::CRS or sf::crs classes?

[rouault]

OGRSpatialReference::importFromWkt() leverages proj_create_from_wkt() which has options to get warnings and errors.

[hobu]

I believe it was Merrick MARS that had a broken WKT writer for a release or two. Is that software_id in the file?

In cases such as these, with PDAL were throwing an error to the user and telling them to manually override them to something resembling what was provided that's actually valid.

[Jean-Romain]

@hobu the header contains

System identifier: NIIRS10
Generating software: GeoCue LAS Updater

[hobu]

@Jean-Romain hmph. I haven't seen bunk WKT out of that software before. Maybe they allowed the user to manually copy/paste WKT into the file or something (RIEGL software allows this!).

[Jean-Romain]

@hobu FYI I retrieved the other file with the WKT mentioned by @rsbivand

System identifier: ALSXX 
Generating software: CloudPro1.2.4.106713

[rsbivand]

I feel that rlas and lidr will need to check WKT credibility (like declared EPSG) - I'll try to provide and expose a function in rgdal to do this, probably based on proj_create_from_wkt().

[rouault]

Slightly related to that topic : #2024 . This will recognize "official" WKT1:ESRI syntax for VERTCS[]. Will not directly help for the 2 cases mentioned in #2022 (comment)

[rsbivand]

Thanks, sounds sensible to accommodate a potentially larger user group.