3

As part of a geo-indexing script, I'd like to add in support for numerous .laz files we have sitting around on our SAN. The only intent of the handler is to grab the filename/path, the SRS, and the extent of the point cloud file. I do not need to do any processing or analysis on the dots.

So I'm able to do

from laspy.file import File

my_file = File("path/to/file.laz", mode='r')

and it loads the files without error. I can call the File.header property, but have no idea how to "get into" it to extract values or calculate the min/max X/Y.

2 Answers 2

5

You should use min and max property of header.

from laspy.file import File

f = File("path/to/file.laz", mode='r')
h = f.header

# h.min: [min_x, min_y, min_z] - h.max: [max_x, max_y, max_z]
extent = [*h.min, *h.max] # extent: [min_x, min_y, min_z, max_x, max_y, max_z]
6
  • Thanks Kadir; is there a way for me to discern the spatial reference from the header? Or just dump all the header contents as a dict and look through them? I was able to grab the min/max as you showed but it was just line and sample, not coordinate information.
    – auslander
    May 22, 2020 at 20:49
  • I've just looked at laspy docs. It says "srs is not implemented". So, you should probably dump the header contents or use liblas. May 22, 2020 at 21:20
  • Great. How do I dump the header contents? That was part of the question alluded to in my original question on how to get into the header to examine all the values.
    – auslander
    May 23, 2020 at 14:54
  • 2
    Use PDAL to get your SRS info. May 29, 2020 at 1:09
  • 2
    @HowardButler how do I use PDAL to get my SRS info?
    – auslander
    Oct 16, 2020 at 16:44
2

As @HowardButler says, you can also use PDAL for this. The documentation is a little hard to find, but you can do it without any need to call subprocess.

First, we need to load the file as a PDAL pipeline:

import pdal
from shapely.geometry import shape, box
from shapely.ops import transform
from pyproj import Transformer
import fiona

filename = "/path/to/las"
pipeline = pdal.Reader.las(filename=data, count=1).pipeline()
pipeline.execute()

count=1 here makes PDAL load only a single datapoint for speed. I would assume it's possible that some LAS files don't contain this metadata in the header and you need to iterate over all points to find the min/max. For my test dataset it doesn't make a difference, but to be safe, omit count.

Then you can access the metadata property to load the info needed (all returned as dictionaries). I'd recommend exploring the dictionary keys here to figure out what info you need, but for example:

# Get the metadata
meta = pipeline.metadata['metadata']

# Bounds
minx = meta['readers.las']['minx']
miny = meta['readers.las']['miny']
maxx = meta['readers.las']['maxx']
maxy = meta['readers.las']['maxy']
las_bbox = box(minx, miny, maxx, maxy)

import pyproj
las_crs = pyproj.crs.CRS(meta['readers.las']['srs']['proj4'])

From here you can do whatever you need - for example, load an arbitrary shapefile and check if any feature intersects the LAS file:

shapefile_path = "/path/to/shapefile"

with fiona.open(shapefile_path, "r") as shapefile:

    project = Transformer.from_crs(shapefile.crs, las_crs, always_xy=True).transform

    for i, feature in enumerate(shapefile):
        geometry = feature["geometry"]
        feature_bounds = transform(project, shape(geometry))
        print(feature_bounds.within(las_bbox))

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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