2

I have older code to crop a chunk within a function (to avoid edge effects when later merging CHM tiles together (i.e., edge pixels derived without enough data. The extend of 1 is to avoid little gaps)), similar to the idea in example 1 here: https://www.rdocumentation.org/packages/lidR/versions/2.2.1/topics/catalog_apply

Raster works but I would like to use Terra instead. I get this error:

Error: Cannot get an Extent object from argument y

I could try to convert the SpatExtent to an Extent, but it would be clunky. Maybe there's an easy way?

Below is part of my code (the Raster lines which work are commented out):

library(raster)
library(terra)

chunk <- ??? 
chm <- ??? 

bbox <- terra::ext(chunk)
bbox <- terra::extend(bbox, 1)
chm  <- terra::crop(chm, bbox)

#bbox <- raster::extent(chunk)
#bbox <- raster::extend(bbox, 1)
#chm  <- raster::crop(chm, bbox)

edits: minimal example, but with lidR code...

library(terra)
library(lidR)

ctgNorm <- readLAScatalog(system.file("extdata", "MixedConifer.laz", package="lidR"))

ctgNorm@output_options$drivers$Raster$param$overwrite <- TRUE

makeCHM_function = function(chunk) {
  las = readLAS(chunk)
  if (is.empty(las))
    return(NULL)
  chm <- grid_canopy(las, 1, pitfree(thresholds = c(0,2,5,10,15), max_edge = c(4, 2.5), subcircle = 0.3))
  rast(chm)
  # Wrong way to get extent... bbox <- terra::ext(chunk)
  # correct way:
  bbox <- lidR::extent(chunk)
  bbox <- terra::extend(bbox, 1)
  chm  <- terra::crop(chm, bbox)
  return(chm)
}

opt_chunk_size(ctgNorm) <- 25
opt_progress(ctgNorm) <- TRUE
opt_output_files(ctgNorm) <- "D:/processingTest20230509a/chm/chmNorm_tile_{XLEFT}_{YBOTTOM}"
myCHM <- catalog_apply(ctgNorm, makeCHM_function)
3
  • 1
    Can you please edit your question and include a minimal, self-contained, reproducible, example that shows the output (error message) where it occurs? That is, add some example data generated by code. I have remove code that you do not need, and added some that you do. May 10, 2023 at 16:59
  • thanks Robert. I had trouble with the sample laz file. I put lidR specific code in if that helps.
    – Ray J
    May 10, 2023 at 17:50
  • In your lidR example, it looks like you are using "raster" code (the use of terra:: is meaningless as explained in my answer). bbox <- terra::ext(chunk) is the way to go if you want to use terra. May 10, 2023 at 18:08

2 Answers 2

2

Your code can simply be written

library(terra)
library(lidR)

ctgNorm <- readLAScatalog(system.file("extdata", "MixedConifer.laz", package="lidR"))

ctgNorm@output_options$drivers$Raster$param$overwrite <- TRUE
opt_chunk_size(ctgNorm) <- 25
opt_progress(ctgNorm) <- TRUE
opt_output_files(ctgNorm) <- "D:/processingTest20230509a/chm/chmNorm_tile_{XLEFT}_{YBOTTOM}"
chm = rasterize_canopy(ctgNorm, 1, pitfree(thresholds = c(0,2,5,10,15), max_edge = c(4, 2.5), subcircle = 0.3))
3
  • Thanks JRR that's certainly better and I'll use that instead of my code. In the case of the tree_detection and rumple_index_surface examples from the link above, would lidR::extent() replace raster::extent()? Or is this better answered in a new question? rdocumentation.org/packages/lidR/versions/2.2.1/topics/…
    – Ray J
    May 10, 2023 at 19:35
  • Ask a new question. You are using old functions deprecated so the best is to start on a new clean basis.
    – JRR
    May 10, 2023 at 19:53
  • 2
    There is no such thing as lidR::extent because extent is a generic function. lidR implements the method for LAS* objects. raster implements it for Raster* objects. These are complimentary implementations of the same method (name) for different objects. One does not replace the other. May 11, 2023 at 1:56
2

I am guessing that chunk and chm are RasterLayer or other Raster* objects.

library(raster)
chm   <- raster(ncol=3, nrow=3, ext=extent(0,2,0,2))
chunk <- raster(ncol=3, nrow=3, ext=extent(0,1,0,1))

Now you do

library(terra)
bbox <- ext(chunk)
bbox <- extend(bbox, 1)
chm  <- crop(chm, bbox)
# Error in .local(x, y, ...) : Cannot get an Extent object from argument y

This happens because "bbox" is a "SpatExtent" and "chm" is a "RasterLayer". The crop method for RasterLayer (implemented in the "raster" package) does not know about SpatExtent, it needs an Extent. To use a SpatExtent, you need to use the method for SpatRaster, implemented in the "terra" package:

chm <- rast(chm)
chm  <- crop(chm, bbox)

It seems that you (with many others) believe that there is a difference between a call to terra::crop(x, y) and to raster::crop(x, y). That is not the case. crop is a generic function, the implementation you get depends on the class of x, not on the namespace::.

Here is an illustration. Below I use the ncol method implemented in "raster" for a RasterLayer --- using base:: does not change anything.

base::ncol(chunk)
#[1] 3

This situation can be different when using standard (non-generic) functions with the same name from different packages, especially when there are name conflicts (including from generics defined independently in different packages).

Bottom line: if you want to use terra methods, and you get Raster* objects from another package, you should start by coercing these to a SpatRaster.

6
  • I think I see that now. bbox isn't a raster. I have to look at this a bit to see why Raster was able to get the extent of bbox...
    – Ray J
    May 10, 2023 at 17:51
  • Okay, I messed that up. Although bbox <- raster::extent(chunk) works, it really should be bbox <- lidR::extent(chunk). I didn't realize lidR had the extent() function too. Sorry about that. I did edit the code above so it works with Terra and Raster isn't needed. Thanks.
    – Ray J
    May 10, 2023 at 18:07
  • That is not correct. lidR::extent returns an Extent; it is the same as raster::extent (see my answer...). If you do that you are using "raster" methods, not "terra" methods. May 10, 2023 at 18:11
  • My question should have been worded better in hindsight. I wanted to run the lidR catalog_apply() without the Raster package in the function. In the old lidR 2.2.1 example, raster::extent() is used. I don't know why lidR::extent() works but I can now just rely on lidR for extent and Terra for extend and crop.
    – Ray J
    May 10, 2023 at 18:26
  • 1
    I just read it again and now I see that my lack of understanding of namespace and dependencies was the root of my confusion. You answered a question I didn't know I had, thanks.
    – Ray J
    May 11, 2023 at 12:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.