2

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.)

3

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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