3

I want to make a simple clip with my .las 3D-Pointcloud.

Therefor I use the lasclip() function.

Import .las:

lidar <- readLAS("~/data.las")
lidar

gives me

class        : LAS (LASF v1.2)
memory       : 774.7 Mb 
extent       : 3403298 , 3403640 , 5286553 , 5286831 (xmin, xmax, ymin, ymax)
area         : 62291.4 m² (convex hull)
points       : 11946448 points,  4704414 pulses
density      : 191.78 points/m²,  75.52 pulses/m²
field names  : X Y Z gpstime Intensity ReturnNumber NumberOfReturns Classification ScanAngle R G B pulseID 
coord. ref.  : +init=epsg:31467 +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +datum=potsdam +units=m +no_defs +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7 

Import the clip:

site <- readOGR(dsn = "~/directory", layer = "clipshape")
site

gives me:

class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : 3403372, 3403590, 5286594, 5286817  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +datum=potsdam +units=m +no_defs +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7 
variables   : 1
names       : id 
value       :  1 

Ok now I want to clip it:

lasclip(lidar, site)

and I get the following Error:

Error: Geometry not supported

I am working in an environment where I cannot use pdal with python for several reasons.

How can I get cliplas() working?

3

Edit: The package is evolving. As stated in JRR's answer which is one of lidR's authors, newer versions (> 1.5) will support clipping with SpatialPolygonDataFrame. It is also planned to support clipping with multipart polygons.


The package ‘lidR’ (version 1.4.2*) pdf says on lasclip:

Usage:
lasclip(x, geometry, ofile = "")

Arguments
geometry: a geometric object. Currently* Polygon from sp is supported.

*apr/2018

Therefore, one needs to clip the .las/.laz data with an object of class Polygon and not SpatialPolygonsDataFrame.

This is a reproducible example which works:

library(lidR)
file.path <- system.file("extdata", "Megaplot.laz", package="lidR")
lidar <- readLAS(file.path)

library(rgdal)
site_spdf <- readOGR(dsn = "...lidR\\extdata", layer = "lake_polygons_UTM17")
site_p <- site_spdf@polygons[[1]]@Polygons[[1]] #sp object of class Polygon

clipped_las = lasclip(lidar, site_p)

Other software options for clipping LiDAR files with vector data are available in Clipping LAS data using shapefile polygons and open source software?

  • Thanks for the answer. I also read the documentation but I could not get the Polygon. site_spdf@polygons[[1]]@Polygons[[1]] is exactly what I was looking for. – Hans Jürgen Apr 23 '18 at 17:17
  • @HansJürgen, no problem. You are welcome. – Andre Silva Apr 23 '18 at 17:39
  • I will accept it as soon as I tried it out, which I cannot do at the moment because the data is not on my private computer. Tomorrow morning I will see if it really solved my problem. – Hans Jürgen Apr 23 '18 at 17:41
  • @HansJürgen, calling the slot object is not necessary, you can coerce on the fly. As pointed out in a comment to your original post you can simple use: as(site_spdf, "Polygons") as the object you pass to the function. – Jeffrey Evans Apr 23 '18 at 17:54
  • 1
    Sorry, in looking at the source for lasclip it is forcing an object class of "Polygon". This is just lazy coding as what is really being called, and passed to lasclipPolygon, is just the coordinates of the polygon. As long as it represents a single object geometry, this could be done with any polygon object class. You can pass the object using lasclip(lidar, Polygon(coordinates(site_spdf[1,]) The reason I do not care for calling the polygons slot directly is that, because the list object will be nested, it will fail if using a multipart sp object. – Jeffrey Evans Apr 23 '18 at 18:22
2

Perfect answer by @andre-silva but let met add few informations. In lidR 1.5.0 you will be able to clip using a SpatialPolygonDataFrame. In that case you will get a list of LAS objects (one per polygon).

las = readLAS("file.las")
spdf <- readOGR(dsn = "...", layer = "...")
clipped_las = lasclip(ctg, spdf)

Also lasclip will be compatible with a LAScatalog. You don't need anymore to load the entire tile to extract a single polygon. Just extract the polygon.

ctg = catalog("directory/")
spdf <- readOGR(dsn = "...", layer = "...")
poly <- spdf@polygons[[1]]@Polygons[[1]] # sp object of class Polygon
clipped_las = lasclip(ctg, poly)

But be careful. You can't clip a SpatialPolygonDataFrame from a LAScatalog. It is technically possible but it has been disabled for two reasons. The first one is memory safety. Indeed, assuming you have 1000 tiles over 1000 km² and a shapefile with hundreds of large polygons, this would load too much data. In order to prevent R crash with bad usage this option is disabled. The second reason is multipart polygons

Regarding multipart polygons, the rlas package only supports single polygon extraction. Thus it is extremely easy to efficiently extract a SpatialPolygon. For SpatialPolygons (potentially multipart polygons with holes) this is not currently natively supported. It requires more computation that are not memory efficient. This is why it is not supported yet but this will be added in lidR later.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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