Error in rasterize_terrain due to projection check when argument res provided as raster

[jmmonnet]

Hi,

I get an error from the rasterize_terrain function when I provide the res argument as a raster and this raster has a projection information. I am not sure I provide the projection information of the raster correctly, see example below. In case the las data and raster have no projection information, then the functions runs correctly.

# install lidaRtRee 3.1.1 package to acces to sample data
# devtools::install_git("https://gitlab.irstea.fr/jean-matthieu.monnet/lidaRtRee/")

library(lidR)
LASfile <- system.file("extdata", "las_chablais3.laz", package="lidaRtRee")
las_chablais3 <- lidR::readLAS(LASfile)

# without projection information
dtm <- raster::raster(xmn = 974330, xmx = 974400, ymn = 6581620, ymx = 6581700, resolution = 1)
dtm <- lidR::rasterize_terrain(las_chablais3, dtm, lidR::tin())

# set projection information
lidR::projection(las_chablais3) <- 2154
dtm <- raster::raster(xmn = 974330, xmx = 974400, ymn = 6581620, ymx = 6581700, resolution = 1, crs = lidR::projection(las_chablais3))
dtm <- lidR::rasterize_terrain(las_chablais3, dtm, lidR::tin())
#> Error in st_crop.stars(x, i, crop = crop, ...): for cropping, the CRS of both objects have to be identical

[Jean-Romain]

Thank you for reporting. This is because the CRS representations between raster, terra and stars are very different and the conversion is sometime not exactly similar. I'll investigate.

[Jean-Romain]

I fixed it. I'm not sure if it is possible to encounter the problem in other functions but it is likely. In short when you are doing:

dtm <- raster::raster(xmn = 974330, xmx = 974400, ymn = 6581620, ymx = 6581700, resolution = 1, crs = lidR::projection(las_chablais3))

The crs of las_chablais3 and dtm are not the same. The first one is a crs with a WKT string. The second is a CRS with a proj4string. And the proj4string is not convertible back to the WKT string from which it originates. The R spatial ecosytem is moved toward WKT CRS this is one of the reasons why I made v4 for lidR and why we should move to sf/stars/terra.

Below some examples that now work

library(lidR)

LASfile <- system.file("extdata", "las_chablais3.laz", package="lidaRtRee")
las <- readLAS(LASfile)
lidR::projection(las) <- 2154
las
#> class        : LAS (v1.2 format 1)
#> memory       : 7 Mb 
#> extent       : 974326, 974408, 6581619, 6581702 (xmin, xmax, ymin, ymax)
#> coord. ref.  : RGF93 / Lambert-93 
#> area         : 6888 m²
#> points       : 92.1 thousand points
#> density      : 13.37 points/m²

library(raster)
rr <- raster(xmn = 974330, xmx = 974400, ymn = 6581620, ymx = 6581700, resolution = 1, crs = lidR::projection(las))
rr
#> class      : RasterLayer 
#> dimensions : 80, 70, 5600  (nrow, ncol, ncell)
#> resolution : 1, 1  (x, y)
#> extent     : 974330, 974400, 6581620, 6581700  (xmin, xmax, ymin, ymax)
#> crs        : +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
dtm <- rasterize_terrain(las, rr, tin())
dtm
#> class      : RasterLayer 
#> dimensions : 80, 70, 5600  (nrow, ncol, ncell)
#> resolution : 1, 1  (x, y)
#> extent     : 974330, 974400, 6581620, 6581700  (xmin, xmax, ymin, ymax)
#> crs        : +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs 
#> source     : memory
#> names      : Z 
#> values     : 1349.18, 1378.4  (min, max)
dtm <- grid_terrain(las, rr, tin())
dtm
#> class      : RasterLayer 
#> dimensions : 80, 70, 5600  (nrow, ncol, ncell)
#> resolution : 1, 1  (x, y)
#> extent     : 974330, 974400, 6581620, 6581700  (xmin, xmax, ymin, ymax)
#> crs        : +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs 
#> source     : memory
#> names      : Z 
#> values     : 1349.18, 1378.4  (min, max)

library(terra)
tr <- rast(rr)
dtm <- rasterize_terrain(las, tr, tin())
dtm
#> class       : SpatRaster 
#> dimensions  : 80, 70, 1  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : 974330, 974400, 6581620, 6581700  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs 
#> source      : memory 
#> name        :       Z 
#> min value   : 1349.18 
#> max value   :  1378.4
terra::crs(tr) <- st_crs(las)$wkt
dtm <- rasterize_terrain(las, tr, tin())
dtm
#> class       : SpatRaster 
#> dimensions  : 80, 70, 1  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : 974330, 974400, 6581620, 6581700  (xmin, xmax, ymin, ymax)
#> coord. ref. : RGF93 / Lambert-93 (EPSG:2154) 
#> source      : memory 
#> name        :       Z 
#> min value   : 1349.18 
#> max value   :  1378.4

library(stars)
sr <- st_as_stars(rr)
dtm <- rasterize_terrain(las, sr, tin())
dtm
#> stars object with 2 dimensions and 1 attribute
#> attribute(s):
#>       Min. 1st Qu.  Median     Mean 3rd Qu.   Max.
#> Z  1349.18 1360.85 1367.98 1366.979 1373.64 1378.4
#> dimension(s):
#>   from to  offset delta                       refsys point values x/y
#> x    1 70  974330     1 +proj=lcc +lat_0=46.5 +lo...    NA   NULL [x]
#> y    1 80 6581700    -1 +proj=lcc +lat_0=46.5 +lo...    NA   NULL [y]
st_crs(sr) <- st_crs(las)
dtm <- rasterize_terrain(las, sr, tin())
dtm
#> stars object with 2 dimensions and 1 attribute
#> attribute(s):
#>       Min. 1st Qu.  Median     Mean 3rd Qu.   Max.
#> Z  1349.18 1360.85 1367.98 1366.979 1373.64 1378.4
#> dimension(s):
#>   from to  offset delta             refsys point values x/y
#> x    1 70  974330     1 RGF93 / Lambert-93    NA   NULL [x]
#> y    1 80 6581700    -1 RGF93 / Lambert-93    NA   NULL [y]