create_split_merge_pdal_crop_and_clip_command.py

#!/usr/bin/python3
"""
Create a shell command to crop point clouds into tiles and perform
classification and clipping by polygon mask on them.

PDAL performs a naive point-in-polygon test for every point in a cloud.
With tens or hundreds of millions of points in a cloud generated by an ODM
run, and a large number of mask polygons, this is unfeasible.

The proper solution is for PDAL itself to crop both the point cloud and 
the mask polygon layer into smaller tiles, perform the point-in-polygon tests
on the individual blocks, and re-merge the resulting classified/tiled point
clouds. While we're at it, we might as well add the option to multi-thread
the operation.

However, not having the time to implement that, here's a workaround hack.

- Create a polygon grid in QGIS. Crop it to your point cloud extent.
- Intersect your mask polygon layer with the grid. Retain the grid coords 
and make sure there's a unique ID.
- Export the mask polygon layer as a CSV, including the tile coords.
- Split the mask polygon layer into multiple files in a folder, using the
grid unique ID as the split coordinate.
- Place this script in a directory containing:
  - A subdirectory called split_mask with the multiple mask files
  - A subdirectory called cropped
  - A subdirectory called clipped
  - The original point cloud file
- Run this script with the folling arguments:
  - The original point cloud file (.las or .laz)
  - The CSV file you exported from mask polygon layer
  - The name of the shell script you want to create
- Set the resulting shell script to be executable and run it.

If you feel like it, you can modify the resulting shell script to multithread
by adding & to the end of each line and slapping in a wait command every so
often (depending how many you want to run concurrently; I'm doing 4 at a time).
When I get around to it, I'll automate that in this script, unless I (or 
someone else) gets around to adding something like this to the PDAL overlays
filter!

Then you can rasterize the resulting point cloud, which is full of holes, and
use a Fill nodata agorithm to fill it in. You might want to do some cleaning
before then, as if there are little building edges in the point cloud a naive
Fill nodata agorithm spans the holes from the raised edges.
"""
import sys, os
import csv

def create_pdal_commands(lasfile, csvfile):
    print(f'Creating command from {csvfile}')
    tiles = list(csv.reader(open(csvfile)))
    header = tiles.pop(0)
    print(f'{csvfile} contains {len(tiles)} chunks')
    commands = []
    classifiedlasfile = lasfile

    if True: # TODO: flag if you want to pre-classify your point cloud
        (laspath, lasext) = os.path.splitext(lasfile)
        classifiedlasfile = f'{laspath}_classified_and_clipped{lasext}'
        commands.append('echo "Classifying point cloud using Simple Morphological Filter"')
        commands.append(f'pdal translate {lasfile} -o {classifiedlasfile} outlier smrf range --filters.outlier.method="statistical" --filters.outlier.mean_k=8 --filters.outlier.multiplier=3.0 --filters.smrf.ignore="Classification[7:7]" --filters.range.limits="Classification[2:2]" --writers.las.compression=true --verbose 4')
    
    for tile in tiles:
        tilenum = tile[0]
        maskfile = f'split_mask/fid_{tilenum}.gpkg'
        croppedfile = f'cropped/{tilenum}.laz'
        clippedfile = f'clipped/{tilenum}.laz'
        xmin = tile[1]
        ymax = tile[2]
        xmax = tile[3]
        ymin = tile[4]
        commands.append(f'echo Cropping and clipping tile {tilenum}.')
        if os.path.isfile(maskfile):
            commands.append(f'pdal translate {classifiedlasfile} -o {croppedfile} crop --filters.crop.bounds="([{xmin},{xmax}],[{ymin},{ymax}])" --writers.las.compression=true --verbose 4')
            commands.append(f'pdal translate {croppedfile} -o {clippedfile} overlay range --filters.overlay.datasource={maskfile} --filters.overlay.column="CLS" --filters.overlay.dimension="Classification" --filters.range.limits="Classification[2:2]" --verbose 4')
        else:
            commands.append(f'{maskfile} not present. '
                           'Putting cropped result in clipped directory')
            commands.append(f'pdal translate {classifiedlasfile} -o {clippedfile} crop --filters.crop.bounds="([{xmin},{xmax}],[{ymin},{ymax}])" --writers.las.compression=true --verbose 4')

    return commands

def write_bash_script(commands, outfile):
    with open(outfile, 'w') as of:
        of.write('#!/bin/bash')
        of.write('\n\n')
        
        for command in commands:
            of.write(command)
            of.write('\n\n')
        
if __name__ == "__main__":
    commands = create_pdal_commands(sys.argv[1], sys.argv[2])
    write_bash_script(commands, sys.argv[3])