I've looked up pdal documentation on how to clip LAS point cloud files:- https://pdal.io/workshop/exercises/analysis/clipping/clipping.html

Referencing this example, and having a quick look around, it seems like there only is an option for one output file. Please correct me if I'm wrong.

I am after the input LAS files being split into four separate output LAS files (following the geometries of the indextile shapefile).

I was able to achieve this with the LASTools command lasclip.exe, however I've reached the limit of the 'unlicenced' version of lasclip:-

lasclip -i C:/a/test/BY17_1401merged_test.las -split -poly C:/a/cropped_topo50.shp -o newfile.las


In the stdout for the above command;
WARNING: unlicensed. over 3 million points. output slightly distorted.
         tiny xyz noise. points permuted. intensity, gps_time, user_data
         and point_source_ID zeroed.

In the lasinfo report for one of the output LAS files. 
WARNING: real max x larger than header max x by 0.020000
WARNING: real min x smaller than header min x by 0.020000
WARNING: real max y larger than header max y by 0.020000
WARNING: real min y smaller than header min y by 0.020000
WARNING: real min z smaller than header min z by 0.010000

Is there any other way of splitting a LAS file into multiple output files using shapefile geometries? I have a Global Mapper license but cannot find a way to do it with this software. Unfortunately I was unable to get SAGA to work on my laptop and Whitebox is almost unusable (very slow) on my ubuntu computer.

  • Are you trying to avoid the licensing restriction in this software.. I recommend you purchase a license, it's not that expensive and supports the development of future products. I had the same problem trying to extract by irregular polygons which left me two options: export a shapefile for each unique polygon and call lasclip on each shapefile or write my own.. I wrote my own utility in C++ as it seemed the less tedious option. – Michael Stimson Nov 20 '18 at 7:04
  • @MichaelStimson thank you for the heads up. As I'm using these tools for a govt department, it is not so cheap in this situation. Are you aware of any open source options to solve this issue? – Rose Nov 20 '18 at 8:43
  • You could write a python script using laspy and OGR, looping through each record and checking within gdal.org/python/osgeo.ogr.Geometry-class.html#Within and writing to a new las file with laspy if True. – Michael Stimson Nov 21 '18 at 0:13
  • Thank you for the suggestion, for OGR_G_Within would that not take into account points that fall on the geometry boundaries, and cause discrepancies between the number of total output points compared to input? I've found a way to split LAS by Shapefile boundaries using FME, conserving points however this is not an open source option, and very cumbersome. Only feasible for one small LAS file. If someone is interested I can put the workflow in the answer dialog. – Rose Nov 24 '18 at 22:44
  • 1
    Within means within, not on a boundary. In reality it's extremely unlikely for a point to fall on a boundary... coordinates in shapefiles are stored as double but las records are integers divided by a scale factor and offset - not the full double range, without getting too technical it would take a perfect storm situation for a LAS record to be on a line segment. – Michael Stimson Nov 25 '18 at 23:15

The lidR package can do that in R


ctg = catalog("C:/a/test/BY17_1401merged_test.las")
shp = shapefile("C:/a/cropped_topo50.shp")

opt_output_file(ctg) = "C:/a/test/BY17_1401_split_{ID}"

new_ctg <- lasclip(BY17, shp)

This script references your file (or your files if several). It loads your shapefile. It assigns a templated output filename. Here your files will be named BY17_1401_split_1.las, BY17_1401_split_2.las and so on. Then it clips all your polygons by streaming the input file(s) into the output files. To finish it loads back a reference to the new created files.

The templated filename can be more elaborated. For example let assume that your table of attribute has a field names "OWNER" because the polygons describe land ownership your can:

opt_output_file(ctg) = "C:/a/test/BY17_1401_split_{ID}_{OWNER}"

Here your files will be named BY17_1401_split_1_bob.las, BY17_1401_split_2_alice.las and so on.

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