r-lidar / lidR

Airborne LiDAR data manipulation and visualisation for forestry application
https://CRAN.R-project.org/package=lidR
GNU General Public License v3.0
596 stars 131 forks source link

lasclip method for spatialpolygonsdataframe. Iterate through list of polygons. #118

Closed bw4sz closed 6 years ago

bw4sz commented 6 years ago

It would be nice to be able to directly grab a list of point clouds from the output of lastrees object. Let me know if something like this exists, if not I can submit a PR.

I'm working with a toy function similar to this.

segment_ITC<-function(las){

  ws = seq(3,21, 3)
  th = seq(0.1, 2, length.out = length(ws))

  lasground(las, "pmf", ws, th)

  # normalization
  lasnormalize(las, method = "knnidw", k = 10L)

  # compute a canopy image
  chm= grid_canopy (las, res=1, subcircle = 0.1, na.fill = "knnidw", k = 3, p = 2)
  chm = as.raster(chm)
  kernel = matrix(1,3,3)
  chm = raster::focal(chm, w = kernel, fun = mean)
  chm = raster::focal(chm, w = kernel, fun = mean)

    ttops = tree_detection(chm, 5, 2)
    crowns<-lastrees_silva(las, chm, ttops, max_cr_factor = 0.6, exclusion = 0.3,
                   extra = T)

    # display
    tree = lasfilter(las, !is.na(treeID))
    plot(tree, color = "treeID", colorPalette = pastel.colors(100), size = 1)

    contour = rasterToPolygons(crowns, dissolve = TRUE)

    plot(chm, col = height.colors(50))
    plot(contour, add = T)
    return(list(polygons=contour,raster=crowns,tops=ttops))
}
LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
tile = readLAS(LASfile, select = "xyz", filter = "-drop_z_below 0")
col = pastel.colors(200)
silva2016<-segment_ITC(tile,algorithm = "silva2016")

What I was hoping for was to be able to clip the las based on the contents of a SpatialPolygonsDataFrame.

lasclipPolygon(tile,silva2016$polygons)

which would return a list of LAS objects.

but docs suggest I need to decompose this into its coords for a single polygon.

The xpoly coords for a sample poylgon (id ==10) are

silva2016$polygons[silva2016$polygons$id==10,]@polygons[[1]]@Polygons[[1]]@coords[,1]

and the y poly coords are

silva2016$polygons[silva2016$polygons$id==10,]@polygons[[1]]@Polygons[[1]]@coords[,2]

We could just iterate through every polygon and save it and return a list? Thoughts?

First few hours with your package, great work @Jean-Romain. Hopefully I can make myself useful in contrib. I think others would be excited for a direct method. This would be my first PR to this package, so give me a sense for conventions (naming, inputs, regenerate the roxygen?).

Jean-Romain commented 6 years ago

You are the second this week to ask me for this feature. This is not complex to do. I can do it on Monday.

But you missed a point. A generic function must also consider multi-part polygons and polygons with holes. It is a bit more complex than simply iterate through polygons. lasclassify already supports SpatialPolygonDataFrame and can be used as basis for such function.

bw4sz commented 6 years ago

sounds good. I am happy to help in anyway. It might also be worth having a flag to spit out a list, or select multiple features. For example,I might want to see a set of polygons (id %in% c(1,2,3)). I can think of a couple ways to do this. As I said, I just started playing with the package. Thanks for the work.

Jean-Romain commented 6 years ago

Btw looking at your code I can't see any obvious interest of such feature in your specific case. If you want to compute some metrics for each tree you can use tree_metrics. And if you really want to split the point cloud you can use lasfilter and filter for each tree ID.

Anyway this feature might be useful. The difficulty is not to program it for LAS objects but also make it compatible in a seamless way for LAScatalog objects.

Jean-Romain commented 6 years ago

Here a ready to use function to split a LAS object using a SpatialPolygonDataFrame

lassplitspdf = function(las, spdf)
{
  id = lidR:::classify_from_shapefile(las, spdf)
  X = split(las@data, id)
  X = lapply(X, LAS, header = las@header)
  return(X)
}
bw4sz commented 6 years ago

Great. This implementation works well because if someone wanted a group of trees, they could assign the ID field to encompass a set.