Extract Point Cloud bounding box with PDAL in Python?

Question

Data:

I have a point cloud file, let call it "pc_file" (it's not a .las but a .ply file ; https://en.wikipedia.org/wiki/PLY_(file_format) ).

Goal to achieve:

I am searching a fast way to extract the bounding box in the xy plane.

To do that, from the working directory, I can:

A)

1. Run pdal info --boundary pc_file > boundary.json

2. Parse the json file through a python script where I need the bbox:

   import os
   from shapely import wkt
   working_dir = '/absolute/path/to/my/WorkingFolder/'
   boundaryFile = os.path.join(working_dir, 'boundary.json')
   with open(boundaryFile) as f:
    data = json.load(f)
   wkt_boundary1 = data[u'boundary'][u'boundary'].encode('ascii','replace')
   bbox1 = wkt.loads(wkt_boundary1).bounds
   

This is clean, but the main drawback is that I need two different commands; one bash line followed by a run of the python script (of course I can stack them both in a bash file but I will end up with a temporary 'boundary' file written on the disk, which I do not want if possible). That was why I was searching to do that all at once without any disk writings:

B)
1. Run an other python script where I directly call the pdal module ( https://pdal.io/python.html ):

 import os
 import numpy as np
 from io import StringIO 
 import pandas as pd
 import pdal
 from shapely import wkt
 working_dir = '/absolute/path/to/my/WorkingFolder/'
 PCFile = os.path.join(working_dir, 'pc_file')
 json_str = u"""
 {
 "pipeline": [
 """+'"'+PCFile+'"'+""",
 {
 "type": "filters.hexbin",
 "threshold":10
 }
 ]
 }"""
 print("json string: {}".format(json_str))
 pipeline = pdal.Pipeline(json_str)
 pipeline.validate()
 pipeline.loglevel = 8 #really noisy
 count = pipeline.execute()
 arrays = pipeline.arrays
 metadata = pipeline.metadata
 log = pipeline.log
 json_obj = pd.read_json(StringIO(metadata))['metadata']['filters.hexbin']
 wkt_boundary = json_obj[1][u'boundary'].encode('ascii','replace')
 bbox = wkt.loads(wkt_boundary).bounds

But this second solution takes really longer as it computes a close envelope-like boundary ( https://pdal.io/stages/filters.hexbin.html ):

Boundary from the **A)** script Boundary from the **B)** script

The question:

I did not found a way to call a "pdal info" like command directly in python from the pdal module (which would be, in my humble opinion, a clean way to extract point cloud metadata from python).
The pdal module seems to only deal with "pipeline" structures, thus calling "filters".
The pdal "info" command does not exist as a filter ( https://pdal.io/stages/filters.html ).

Is there any or some filter to extract the bounding box directly?

(Don't hesitate to tell me if the question would have better chances on SO directly.)

Answer 1

Seems like you should get to know Python's built-in subprocess module. Among other things, it enables you to make command line calls from inside a script and capture the results.

import subprocess
import json
from shapely.geometry import Polygon
result = subprocess.run(['pdal', 'info', '/path/to/file.ply'],
 stderr = subprocess.PIPE, # stderr and stdout get
 stdout = subprocess.PIPE) # captured as bytestrings
# decode stdout from bytestring and convert to a dictionary
json_result = json.loads(result.stdout.decode())

From here, you can inspect all the information that `pdal info` returned inside your `json_results` dictionary, including getting info about the bounding box.

Here's what json_results['stats']['bbox']['native'] returned for an example lidar file I have:

{'bbox': {'maxx': 552499.99,
 'maxy': 5187999.99,
 'maxz': 349.43,
 'minx': 552250,
 'miny': 5187750,
 'minz': 250.86},
 'boundary': {'coordinates': [[[552250.0, 5187750.0],
 [552250.0, 5187999.99],
 [552499.99, 5187999.99],
 [552499.99, 5187750.0],
 [552250.0, 5187750.0]]],
 'type': 'Polygon'}}

If you want to go straight to making a polygon, you could do something like this:

coords = json_results['stats']['bbox']['native']['boundary']['coordinates']
bbox_poly = Polygon(*coords)

Keep in mind that these coordinates are in whatever projection your data is in. This projection information is not carried along with the polygon, so you may need to add projection info using another package (like `Fiona` or `GeoPandas`) if you want to export this bounding box polygon to a file.