0

I have two rasters that are the result of the grid_density() function from the liDR package. With those rasters I want to compare point distribution of subsets of a LAS file. My problem is, that I get rasters with different extent (since grid_density() aligns the extent of the output raster exactly to the point cloud. If the point cloud is filtered, points at the edges might be missing and the output raster shrinks a bit.

Often, this is only on one side and by a few rows or colums. Since the rasters are generally in the same resolution, I'd rather avoid a resample(), since it is slow and requires interpolation that shouldn't be needed in most cases.

different extents

A code example:

# create reference raster with all points
dgrid_ref = grid_density(las, 2)

# filter las pointcloud and create raster
las_c1 = filter_poi(lasp, Z>=0, Z<5)
dgrid_c1  = grid_density(las_c1, 2)

las_c2 = filter_poi(lasp, Z>=5, Z<9)
dgrid_c2  = grid_density(las_c2, 2)

I tried simply overriding the extent with dgrid_c1 = setExtent(dgrid_c1, dgrid_ref, keepres = TRUE), however that leaves me with a raster that contains only NAs.

I also thought to create an "empty" raster of the reference extent and then "inserting" the subset raster, but couldn't find a way to "paste" a raster into another.

Is there an option to use a "snap raster" for grid_density() that determines the output? Or is there a better option than to resample()? Like overwriting the extent?

2 Answers 2

1

I think you need extend here. Let's make a sample...

Full-world raster:

> rfull = raster()
> rfull[]=1:ncell(rfull)

And extract a smaller raster with the same resolution and origin, so the cells overlap perfectly:

> e = extent(-20,10,-30,0)
> rpart = crop(rfull, e)
> rpart[] = 1:ncell(rpart)

Now we have these two rasters:

> extent(rfull)
class      : Extent 
xmin       : -180 
xmax       : 180 
ymin       : -90 
ymax       : 90 
> extent(rpart)
class      : Extent 
xmin       : -20 
xmax       : 10 
ymin       : -30 
ymax       : 0 

We want to extend the small raster to the extent of the large one, by padding with NA:

> 
> rpack = extend(rpart, rfull)
> 
> extent(rpack)
class      : Extent 
xmin       : -180 
xmax       : 180 
ymin       : -90 
ymax       : 90 

Plotting plot(rpack) gives this:

enter image description here

where you can see the extent is padded out to the extent of the full world raster.

(Note that using extend on a big raster to the extent of a small raster results in an unchanged raster.)

2
  • cool, I didn't know about this function. Thanks a lot!
    – Honeybear
    Commented Sep 7, 2021 at 18:43
  • Me neither. But I thought "hmmm there must be a function to do this in the raster package...", typed ls(package:raster) and guessed...
    – Spacedman
    Commented Sep 7, 2021 at 18:55
0

As I found out here, there is the option to pass a RasterLayer to the grid_density function (res = raster(xy)), which will then be used as the Layout of the output (e.g. see doc on grid_canopy).

Using the same RasterLayer Layout on both

# create reference raster with all points
dgrid_ref = grid_density(las, 2)

# filter las pointcloud and create raster
las_c1 = filter_poi(lasp, Z>=0, Z<5)
dgrid_c1  = grid_density(las_c1, res=dgrid_ref)

las_c2 = filter_poi(lasp, Z>=5, Z<9)
dgrid_c2  = grid_density(las_c2, res=dgrid_ref)

will then use dgrid_ref for extent and resolution. extend won't be needed then.

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.