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.

Ian Turton
  • 81,417
  • 6
  • 84
  • 185
auslander
  • 1,507
  • 1
  • 11
  • 36

2 Answers2

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]
Kadir Şahbaz
  • 76,800
  • 56
  • 247
  • 389
  • 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 '20 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. – Kadir Şahbaz May 22 '20 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 '20 at 14:54
  • Please check these posts: https://gis.stackexchange.com/a/203395/29431 and https://gis.stackexchange.com/a/80762/29431 – Kadir Şahbaz May 23 '20 at 20:48
  • 2
    Use PDAL to get your SRS info. – Howard Butler May 29 '20 at 01:09
  • 2
    @HowardButler how do I use PDAL to get my SRS info? – auslander Oct 16 '20 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))
Josh
  • 151
  • 4