2

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

https://www.cs.unc.edu/~isenburg/lastools/download/lasclip_README.txt

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.

6
  • 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. Commented Nov 20, 2018 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
    Commented Nov 20, 2018 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. Commented Nov 21, 2018 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
    Commented Nov 24, 2018 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. Commented Nov 25, 2018 at 23:15

1 Answer 1

2

The lidR package can do that in R

library(lidR)
library(sf)

ctg = readLAScatalog("C:/a/test/BY17_1401merged_test.las")
shp = st_read("C:/a/cropped_topo50.shp")

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

# with clip_roi
new_ctg <- clip_roi(ctg, 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.

3
  • it sounds cool, but I cannot reproduce this. I get "cannot find function opt_output_file=. Any ideas? Also I'm thinking the provided code isn't right: The line new_ctg <- lasclip(BY17, shp) probably shouldn't contain BY17 but ctg? And isn't the .las extension missing in the filename template?
    – Honeybear
    Commented Sep 8, 2021 at 22:00
  • ask a new question
    – JRR
    Commented Sep 8, 2021 at 22:46
  • I am not sure a new question is justified here? The question I am asking in my comment is specifically about your answer here, which doesn't work as is. Asking a new question about an answer in a question which is already marked as duplicate doesn't seem the right thing to do here. After some research I found that it you code has a typo opt_output_files(ctg) was missing an s. I made an edit, please accept it.
    – Honeybear
    Commented Sep 13, 2021 at 7:38

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