2

I can do this:

import json
import subprocess as sp

laz_file = somefile.laz
r = (sp.run(['pdal', 'info', laz_file], stderr=sp.PIPE, stdout=sp.PIPE))
json_info = json.loads(r.stdout.decode())

and parse through the JSON, get the info I need (EPSG code and min/max x/y). It does "work."

However calling this workflow on a laz file, especially when some of mine are over 1 gb, takes significantly longer than a las file. As I am iterating over hundreds of laz files (which change daily and sometimes hourly) and intending the output to be used as a dynamic map tracker for our holdings.

How can I extract the EPSG code and the bounding box in a more performant way?

4
  • Where's the bottleneck? If you're IO bound, use faster storage if you can. If you're not IO or memory bound, try using multiprocessing. You could also see if lasinfo is any quicker than pdal, if it provides the info you need.
    – mikewatt
    Commented Oct 21, 2020 at 23:05
  • 1
    Have you worked with pgpointcloud? Check out PC_Summary() github.com/pgpointcloud/pointcloud
    – Aaron
    Commented Oct 22, 2020 at 3:43
  • @mikewatt it's the pdal --stats routine itself, but only on laz files. Run against very large las files, it is quick. So it seems like there is some decompress-to-memory going on that is very slow. I am running this all across 10gb fibre to nvme isilon, so I can't imagine storage IO is the bottleneck.
    – auslander
    Commented Oct 22, 2020 at 14:44
  • 1
    @Aaron that is really intriguing but as I am in a tightly controlled enterprise environment that does not run PG, I will need solutions that let me interrogate data at rest on the san without postgresql.
    – auslander
    Commented Oct 22, 2020 at 14:45

1 Answer 1

2

According to the doc of pdal info:

If no options are provided, --stats is assumed.

This means that you need to read the entire point cloud to derive information that you don't actually need. What you are looking for is in the metadata of the file (the header). You don't need to read the payload (the points) at all. Try

pdal info file.las --metadata 

This will return your info (epsg, bbox) virtually instantaneously. You will maybe need a bit more parsing however. The EPSG code is stored is the Variable Length Record (VLR) with a key "LASF_Projection". Should look like that.

"vlr_0":
{
  "data": "AQABAAAABQAABAAAAQABAAAMAAABAIkLBAwAAAEAKSMDEAAAAQApIwAQAAABAGcE",
  "description": "by LAStools of rapidlasso GmbH",
  "record_id": 34735,
  "user_id": "LASF_Projection"
},

The bounding box can be found at the begining

"maxx": 684993.29,
"maxy": 5018007.25,
"maxz": 29.97,
"minor_version": 2,
"minx": 684766.39,
"miny": 5017773.08,
"minz": 0,
9
  • You are partially right. It returns the XY extents, but only as line/sample, not as geocoordinates like --stats does. And it appears to return the spatial reference as a huge text brick inside of one element - the SRS info itself not being serialized. I'm probably missing something.
    – auslander
    Commented Oct 22, 2020 at 14:32
  • Does my edit clarifies the anwser?
    – JRR
    Commented Oct 22, 2020 at 14:45
  • So is the record_id value the EPSG code?
    – auslander
    Commented Oct 22, 2020 at 14:49
  • Yes. For LAS file < 1.4 at least. For LAS file 1.4 the CRS is recorded as WKT string.
    – JRR
    Commented Oct 22, 2020 at 14:52
  • Ah...this looks really promising, thank you! I'll run a few tests and report back.
    – auslander
    Commented Oct 22, 2020 at 14:55

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.