R lidR does not read the CRS of my las file

Question

I have a LAS 1.4 file that I can read in lidR but the CRS is not read:

las = readLAS("LD31201353.las")
las
#> class        : LAS (v1.4 format 6)
#> memory       : 172.5 Mb 
#> extent       : 3121464, 3123000, 1353000, 1356000 (xmin, xmax, ymin, ymax)
#> coord. ref.  : NA 
#> area         : 4.5 kunits²
#> points       : 1.97 million points
#> density      : 0.44 points/units²

If I look at the header the VLR is empty

las@header
#> File signature:           LASF 
#> File source ID:           0 
#> Global encoding:
#>  - GPS Time Type: Standard GPS Time 
#>  - Synthetic Return Numbers: no 
#>  - Well Know Text: CRS is WKT 
#>  - Aggregate Model: false 
#> Project ID - GUID:        00000001-0012-07e3-0000-000000000000 
#> Version:                  1.4
#> System identifier:        Merrick 
#> Generating software:      MARS 
#> File creation d/y:        188/2018
#> header size:              375 
#> Offset to point data:     377 
#> Num. var. length record:  0 
#> Point data format:        6 
#> Point data record length: 30 
#> Num. of point records:    1965664 
#> Num. of points by return: 1634171 285040 42901 3469 82 1 0 0 0 0 0 0 0 0 0 
#> Scale factor X Y Z:       0.001 0.001 0.001 
#> Offset X Y Z:             3122231 1354500 11639 
#> min X Y Z:                3121464 1353000 11008.84 
#> max X Y Z:                3123000 1356000 11782.48 
#> Variable length records:  void

But the file actually contains a CRS because lasinfo from LAStools can read it

extended variable length header record 1 of 2:
 reserved             0
 user ID              'LASF_Projection'
 record ID            2112
 length after header  1213
 description          'OGC COORDINATE SYSTEM WKT'
   OGC COORDINATE SYSTEM WKT:
   COMPD_CS["NAD83(2011) / Colorado Central (ftUS) + NAVD88 height (ftUS) - US Geoid Model of 2012 B",PROJCS["NAD83(2011) / Colorado Central (ftUS)",PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["central_meridian",-105.5],PARAMETER["latitude_of_origin",37.833333333333],PARAMETER["standard_parallel_1",39.75],PARAMETER["standard_parallel_2",38.45],PARAMETER["false_easting",3000000],PARAMETER["false_northing",1000000],UNIT["US Survey Foot",0.304800609601222,AUTHORITY["EPSG","9003"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","6428"]],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["US survey foot",0.3048006096012192, AUTHORITY["EPSG","9003"]],AXIS["Geodetic longitude",EAST],AXIS["Geodetic latitude",NORTH],AUTHORITY["EPSG","6318"]],VERT_CS["NAVD88 height (ftUS) - US Geoid Model of 2012 B",VERT_DATUM["North American Vertical Datum 1988",AUTHORITY["EPSG","5103"]],UNIT["US Survey Foot",0.304800609601222,AUTHORITY["EPSG","9003"]],AUTHORITY["EPSG","6360"]],PARAMETER[Vertical_Shift,0.0],PARAMETER[Direction,1.0]]

What do I need to do in R to expose the CRS (besides calling lastools from R...)?

Answer 1

The package rlas that is used to read and write las/laz files has a partial support of file format 1.4. In your case the WKT string is stored in the extended variable length record

extended variable length header record 1 of 2:
  reserved             0
  user ID              'LASF_Projection'
  record ID            2112
  length after header  1213
  description          'OGC COORDINATE SYSTEM WKT'
    OGC COORDINATE SYSTEM WKT:
    COMPD_CS["NAD83(2011) / Colorado Central (ftUS) + NAVD88 height (ftUS) - US Geoid Model of 2012 B",PROJCS["NAD83(2011) / Colorado Central (ftUS)",PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["central_meridian",-105.5],PARAMETER["latitude_of_origin",37.833333333333],PARAMETER["standard_parallel_1",39.75],PARAMETER["standard_parallel_2",38.45],PARAMETER["false_easting",3000000],PARAMETER["false_northing",1000000],UNIT["US Survey Foot",0.304800609601222,AUTHORITY["EPSG","9003"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","6428"]],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["US survey foot",0.3048006096012192, AUTHORITY["EPSG","9003"]],AXIS["Geodetic longitude",EAST],AXIS["Geodetic latitude",NORTH],AUTHORITY["EPSG","6318"]],VERT_CS["NAVD88 height (ftUS) - US Geoid Model of 2012 B",VERT_DATUM["North American Vertical Datum 1988",AUTHORITY["EPSG","5103"]],UNIT["US Survey Foot",0.304800609601222,AUTHORITY["EPSG","9003"]],AUTHORITY["EPSG","6360"]],PARAMETER[Vertical_Shift,0.0],PARAMETER[Direction,1.0]]

This part of the header is not read by rlas. This is because I've never encountered such file yet. Now I've an example file I will work on the support of this format. But for the moment the only solution is to manually set the CRS. Yet the WKT string does not look valid. The following throw an error:

sf::st_crs('COMPD_CS["NAD83(2011) / Colorado Central (ftUS) + NAVD88 height (ftUS) - US Geoid Model of 2012 B",PROJCS["NAD83(2011) / Colorado Central (ftUS)",PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["central_meridian",-105.5],PARAMETER["latitude_of_origin",37.833333333333],PARAMETER["standard_parallel_1",39.75],PARAMETER["standard_parallel_2",38.45],PARAMETER["false_easting",3000000],PARAMETER["false_northing",1000000],UNIT["US Survey Foot",0.304800609601222,AUTHORITY["EPSG","9003"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","6428"]],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["US survey foot",0.3048006096012192, AUTHORITY["EPSG","9003"]],AXIS["Geodetic longitude",EAST],AXIS["Geodetic latitude",NORTH],AUTHORITY["EPSG","6318"]],VERT_CS["NAVD88 height (ftUS) - US Geoid Model of 2012 B",VERT_DATUM["North American Vertical Datum 1988",AUTHORITY["EPSG","5103"]],UNIT["US Survey Foot",0.304800609601222,AUTHORITY["EPSG","9003"]],AUTHORITY["EPSG","6360"]],PARAMETER[Vertical_Shift,0.0],PARAMETER[Direction,1.0]]')

You can maybe use the epsg code 2232 according to your other question but I don't know if it is strictly the same than the one in the header. They look pretty similar but I'm definitively not knowledgeable about CRS.

projection(las) <- "EPSG:2232"
sf::st_crs(las)
#> Coordinate Reference System:
#>   User input: NAD83 / Colorado Central (ftUS) 
#>   wkt:
#> PROJCRS["NAD83 / Colorado Central (ftUS)",
#>     BASEGEOGCRS["NAD83",
#>         DATUM["North American Datum 1983",
#>             ELLIPSOID["GRS 1980",6378137,298.257222101,
#>                 LENGTHUNIT["metre",1]]],
#>         PRIMEM["Greenwich",0,
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>         ID["EPSG",4269]],
#>     CONVERSION["SPCS83 Colorado Central zone (US Survey feet)",
#>         METHOD["Lambert Conic Conformal (2SP)",
#>             ID["EPSG",9802]],
#>         PARAMETER["Latitude of false origin",37.8333333333333,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8821]],
#>         PARAMETER["Longitude of false origin",-105.5,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8822]],
#>         PARAMETER["Latitude of 1st standard parallel",39.75,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8823]],
#>         PARAMETER["Latitude of 2nd standard parallel",38.45,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8824]],
#>         PARAMETER["Easting at false origin",3000000,
#>             LENGTHUNIT["US survey foot",0.304800609601219],
#>             ID["EPSG",8826]],
#>         PARAMETER["Northing at false origin",1000000,
#>             LENGTHUNIT["US survey foot",0.304800609601219],
#>             ID["EPSG",8827]]],
#>     CS[Cartesian,2],
#>         AXIS["easting (X)",east,
#>             ORDER[1],
#>             LENGTHUNIT["US survey foot",0.304800609601219]],
#>         AXIS["northing (Y)",north,
#>             ORDER[2],
#>             LENGTHUNIT["US survey foot",0.304800609601219]],
#>     USAGE[
#>         SCOPE["Engineering survey, topographic mapping."],
#>         AREA["United States (USA) - Colorado - counties Arapahoe; Chaffee; Cheyenne; Clear Creek; Delta; Denver; Douglas; Eagle; El Paso; Elbert; Fremont; Garfield; Gunnison; Jefferson; Kit Carson; Lake; Lincoln; Mesa; Park; Pitkin; Summit; Teller."],
#>         BBOX[38.14,-109.06,40.09,-102.04]],
#>     ID["EPSG",2232]]

Edit: with rlas >= 1.3.10 and lidR >= 3.1.2 the extended variable length records are now read and understood by R