point-in-las

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Author: Daniel Rode
# Dependencies:
#   Python 3.11+
#   PDAL
#   shapely
# Created: 28 May 2024
# Updated: -
# Version: 0


# Description: Given a set of coordinates, determine whether the point falls
# within a given las/laz file point cloud boundary.


import sys
from sys import exit

import subprocess as sp
import json

import shapely.wkt
from shapely.geometry import Point


# Constants
exe_name = sys.argv[0].split('/')[-1]  # This script's filename
help_text = f"Usage: python3 {exe_name} X_COORD Y_COORD LAS_PATH..."


# Functions
def get_las_bound(las_path):
    cmd = ['pdal', 'info', '--boundary', las_path]
    p = sp.run(cmd, check=True, capture_output=True, text=True)

    return json.loads(p.stdout)['boundary']['boundary']
    

# Main
def main():
    # Parse command line arguments
    args = sys.argv[1:]
    try:
        x_coord = args[0]
        y_coord = args[1]
    except IndexError:
        print(help_text)
        exit(1)
        
    # Check which las/laz file contains the given coordinate pair
    count = 0
    for las_path in args[2:]:
        bound = shapely.wkt.loads(get_las_bound(las_path))
        point = Point(x_coord, y_coord)
        if bound.contains(point):
            print(las)
            exit()

        print(f"checked {count}")
            

if __name__ == '__main__':
    main()