2

I want to check if the extent of a LAS file overlap any polygon in a shapefile:

library(lidR)
library(rgdal)
library(sf)
library(sp)

f = "D:/Skapa_DEM/inputs/19D021_67000_5250_25.laz"

las <- readLAS(f)
lasextent = extent(las)
xmin = lasextent[1]
xmax = lasextent[2]
ymin = lasextent[3]
ymax = lasextent[4]

e <- as(raster::extent(ymin, ymax, xmin, xmax), "SpatialPolygons")
proj4string(e) <- "+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +type=crs"
df <- readOGR("D:/Skapa_DEM/testarea.shp")

I've tried:

#st_intersection(e, df)
#Error in UseMethod("st_intersection") : 
#  no applicable method for 'st_intersection' applied to an object of class "c('SpatialPolygons', 'Spatial', 'SpatialVector')"

#intersect(df, e)
#Error in as.vector(x) : no method for coercing this S4 class to a vector

But it does not work, see errors. There must be some easier way. How?

1
  • well, your error refers to st_intersection being created for class sf (from package sf) and you are using it on sp objects, you may use sf::st_as_sf() to convert your sp objects and then just try st_intersects (which just returns a logical) of st_intersection if you do need the area extent of the intersection
    – Elio Diaz
    Nov 15, 2022 at 15:36

1 Answer 1

1

You can use functions from the sf package for this.

First, get the LAS extent as an sf polygon (I'm using the output from help(readLAS) here):

> las_extent = st_as_sfc(st_bbox(las))

Then if you have sf polygons with the right CRS:

> pols
Geometry set for 3 features 
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 684733.5 ymin: 5017723 xmax: 684937.2 ymax: 5017947
Projected CRS: NAD83 / UTM zone 17N
POLYGON ((684733.9 5017744, 684738.3 5017772, 6...
POLYGON ((684797 5017926, 684812.9 5017883, 684...
POLYGON ((684881.7 5017937, 684880 5017891, 684...

you can use st_intersects to test intersection:

> st_intersects(pols, las_extent)
Sparse geometry binary predicate list of length 3, where the predicate
was `intersects'
 1: (empty)
 2: 1
 3: 1

This shows my first polygon doesn't intersect the LAS boundary but the other two do.

> plot(las_extent, lwd=2, col="grey")
> plot(pols, add=TRUE)

enter image description here

Grey is my LAS extent, the first polygon is bottom left. The other two intersect.

Doing the test this way you can test the length of each intersection vector using lengths to get a boolean test of intersection:

> lengths(st_intersects(pols, las_extent)) > 0
[1] FALSE  TRUE  TRUE

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.