projection not set on LAScatalog with a crs object from sf

[lucas-johnson]

Current Behavior

When using projection() to set the crs on a LAScatalog the result is NA and no warning/error is produced.

Example:

project_epsg
#> [1] 3625
test <- st_crs(project_epsg)
test
#> Coordinate Reference System:
#>   User input: EPSG:3625 
#>   wkt:
#> PROJCS["NAD83(NSRS2007) / New York East",
#>     GEOGCS["NAD83(NSRS2007)",
#>         DATUM["NAD83_National_Spatial_Reference_System_2007",
#>             SPHEROID["GRS 1980",6378137,298.257222101,
#>                 AUTHORITY["EPSG","7019"]],
#>             TOWGS84[0,0,0,0,0,0,0],
#>             AUTHORITY["EPSG","6759"]],
#>         PRIMEM["Greenwich",0,
#>             AUTHORITY["EPSG","8901"]],
#>         UNIT["degree",0.0174532925199433,
#>             AUTHORITY["EPSG","9122"]],
#>         AUTHORITY["EPSG","4759"]],
#>     PROJECTION["Transverse_Mercator"],
#>     PARAMETER["latitude_of_origin",38.83333333333334],
#>     PARAMETER["central_meridian",-74.5],
#>     PARAMETER["scale_factor",0.9999],
#>     PARAMETER["false_easting",150000],
#>     PARAMETER["false_northing",0],
#>     UNIT["metre",1,
#>         AUTHORITY["EPSG","9001"]],
#>     AXIS["X",EAST],
#>     AXIS["Y",NORTH],
#>     AUTHORITY["EPSG","3625"]]
projection(ctg) <- test
ctg
#> class       : LAScatalog (v1.4 format 6)
#> extent      : 114000, 115500, 580500, 581495.6 (xmin, xmax, ymin, ymax)
#> coord. ref. : NA 
#> area        : 1.49 kunits²
#> points      : 10.72million points
#> density     : 7.2 points/units²
#> num. files  : 1 

Desired behavior

projection() should correctly set the crs of the LAScatalog, but if the projection() results in NA crs then a warning or error message should be produced.

[Jean-Romain]

So I checked the source code and surprisingly there is no method projection<- for LAScatalog objects. What does it means? It means that a LAScatalog is a SpatialPolygonDataFrame and thus this the is the raster::projection() methods that is called. You can reproduce your issue without lidR

library(sf)
library(raster)
library(sp)

x_coord <- c(16.48438,  17.49512,  24.74609, 22.59277, 16.48438)
y_coord <- c(59.736328125, 55.1220703125, 55.0341796875, 61.142578125, 59.736328125)
xym <- cbind(x_coord, y_coord)
p = Polygon(xym)
ps = Polygons(list(p),1)
sps = SpatialPolygons(list(ps))

sps
#> class       : SpatialPolygons 
#> features    : 1 
#> extent      : 16.48438, 24.74609, 55.03418, 61.14258  (xmin, xmax, ymin, ymax)
#> crs         : NA

projection(sps)  <- st_crs(3625)
sps
#> class       : SpatialPolygons 
#> features    : 1 
#> extent      : 16.48438, 24.74609, 55.03418, 61.14258  (xmin, xmax, ymin, ymax)
#> crs         : NA

test <- CRS(SRS_string = "EPSG:3625")
projection(sps) = test
sps
#> class       : SpatialPolygons 
#> features    : 1 
#> extent      : 16.48438, 24.74609, 55.03418, 61.14258  (xmin, xmax, ymin, ymax)
#> crs         : +proj=tmerc +lat_0=38.8333333333333 +lon_0=-74.5 +k=0.9999 +x_0=150000 +y_0=0 +ellps=GRS80 +units=m +no_defs

^{Created on 2021-01-04 by the reprex package (v0.3.0)}

That being said you can do

projection(ctg) <- 3625
#> Warning message:
#> In showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
#>   Discarded datum NAD83 (National Spatial Reference System 2007) in CRS definition
ctg
#> class       : LAScatalog (v1.2 format 1)
#> extent      : 684766.4, 684993.3, 5017773, 5018007 (xmin, xmax, ymin, ymax)
#> coord. ref. : NAD83(NSRS2007) / New York East 
#> area        : 53133.17 m²
#> points      : 81.6thousand points
#> density     : 1.5 points/m²
#> num. files  : 1 

I keep the issue opened for an enhancement.

[lucas-johnson]

It looks like this problem also exists for LAS objects. And I am seeing some different results than what you are showing with your last solution (projection(ctg) <- 36251). Perhaps this is a result of different spatial libraries (GDAL PROJ etc)? Or maybe I just need the latest development version of lidR? I'm currently using 3.0.5

Different outcomes with LAScatalog - your last proposed solution does not work:

ctg
#> class       : LAScatalog (v1.4 format 6)
#> extent      : 167516.2, 167609.8, 617820.2, 617914.3 (xmin, xmax, ymin, ymax)
#> coord. ref. : NA 
#> area        : 8802.14 units²
#> points      : 26.9thousand points
#> density     : 3.1 points/units²
#> num. files  : 1 
projection(ctg) <- 3625
#> Error in CRS(SRS_string = x) : no arguments in initialization list
packageVersion("lidR")
#> [1] ‘3.0.5’

Problem with LAS object as well, however your solution works in this case:

a_las
#> class        : LAS (v1.4 format 6)
#> memory       : 2.4 Mb 
#> extent       : 167516.2, 167609.8, 617820.2, 617914.3 (xmin, xmax, ymin, ymax)
#> coord. ref.  : NA 
#> area         : 6911.356 units²
#> points       : 26.9 thousand points
#> density      : 3.89 points/units²
projection(a_las) <- st_crs(3625)
#> Error in wkt == "" : 
#>   comparison (1) is possible only for atomic and list types
projection(a_las) <- 3625
a_las
#> class        : LAS (v1.4 format 6)
#> memory       : 2.4 Mb 
#> extent       : 167516.2, 167609.8, 617820.2, 617914.3 (xmin, xmax, ymin, ymax)
#> coord. ref.  : +proj=tmerc +lat_0=38.83333333333334 +lon_0=-74.5 +k=0.9999 +x_0=150000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
#> area         : 6911.356 m²
#> points       : 26.9 thousand points
#> density      : 3.89 points/m²

[Jean-Romain]

Please show

rgdal::getGDALVersionInfo()
#> [1] "GDAL 3.2.0, released 2020/10/26"
rgdal::getPROJ4VersionInfo()
#> [1] "Rel. 7.2.0, November 1st, 2020, [PJ_VERSION: 720]"
packageVersion("rgdal")
#> [1] ‘1.5.18’
packageVersion("sf")
#> [1] ‘0.9.6’

I'm currently using 3.0.5

I released 3.1.0 this morning. It should be available on CRAN within few days. Anyway what we are discussing right now is not related to the upcoming modifications because we are talking about stuff introduced in 3.0.4

[lucas-johnson]

rgdal::getGDALVersionInfo()
#> [1] "GDAL 2.2.3, released 2017/11/20"
rgdal::getPROJ4VersionInfo()
#> [1] "Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]"
#> attr(,"short")
#> [1] 493
packageVersion("rgdal")
#> [1] ‘1.5.18’
packageVersion("sf")
#> [1] ‘0.9.6’

Looks like PROJ and GDAL are out of date on my machine...

[Jean-Romain]

Yes. This is definitively a complex topic on which I don't understand everything myself. You can see this question and more importantly the links that are given in the answer.

From my point of view, lidR has been updated to match with the latest development of the spatial ecosystem because it is a requirement to maintain the package on CRAN. There are some conditional statements in the source code to maintain backward compatibility but the truth is that they are not tested and I probably missed some use cases. On CRAN every R flavour use the recent version of GDAL and PROJ, continuous integration as well and on my machine I won't deinstall and reinstall outdated stuff to perform some check (I did it at the beginning).

Consequently my best advice is to be up-to-date because I can hardly figure out and fix stuff related to GDAL < 3 and PROJ < 6

[lucas-johnson]

Ok thanks! Makes sense.

[Jean-Romain]

I added dedicated functions in 3.2.0

library(lidR)
library(sf)

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
ctg = readLAScatalog(LASfile)
ctg
#> class       : LAScatalog (v1.2 format 1)
#> extent      : 684766.4, 684993.3, 5017773, 5018007 (xmin, xmax, ymin, ymax)
#> coord. ref. : NAD83 / UTM zone 17N 
#> area        : 53133.17 m²
#> points      : 81.6thousand points
#> density     : 1.5 points/m²
#> num. files  : 1
projection(ctg) <- st_crs(3625)
ctg
#> class       : LAScatalog (v1.2 format 1)
#> extent      : 684766.4, 684993.3, 5017773, 5018007 (xmin, xmax, ymin, ymax)
#> coord. ref. : NAD83(NSRS2007) / New York East 
#> area        : 53133.17 m²
#> points      : 81.6thousand points
#> density     : 1.5 points/m²
#> num. files  : 1

^{Created on 2021-05-25 by the reprex package (v2.0.0)}