Add PDAL to comparison matrix

[gadomski]

As a part of openjournals/joss-reviews#7123, I used both startinpy and PDAL to create a terrain mesh from a USGS 3DEP point cloud. PDAL uses delaunator-cpp for its triangulation, wrapped in filters.delaunay. Because PDAL has a Python API and has good performance, it might be good to include it in the comparisons. I can't add it to the comparisons myself because they use @hugoledoux local files. EDIT: I see that there's links to download the performance testing files, so I'm going to take a stab at adding PDAL myself 🙇🏼.

I don't want to assume that I'm doing things correctly, but in my tests I found that PDAL + delaunator-cpp was significantly faster in a simple build-terrain operation. I did my best to ensure that I was using a release profile when building startinpy but of course could have made an error.

cc @kylemann16 who is a PDAL developer.

Comparison

I compared performance using a tiled COPC, fetched from Microsoft's Planetary Computer using stac-cli, jq, and stac-asset:

stac search https://planetarycomputer.microsoft.com/api/stac/v1 \
    -c 3dep-lidar-copc \
    --intersects '{"type":"Point","coordinates":[-105.119,40.173]}' \
    | jq '.features[0]' \
    | stac-asset download

startinpy

Code

import time

import laspy
import numpy

import startinpy

start = time.time()
las = laspy.read("USGS_LPC_CO_SoPlatteRiver_Lot2a_2013_13TDE489446_LAS_2015.copc.laz")
points = las.points[las.classification == 2]
points = numpy.stack((points.x, points.y, points.z)).transpose()

read = time.time()
dt = startinpy.DT()
dt.insert(points)

created = time.time()
dt.write_ply("dt-startinpy.ply")

done = time.time()

print(f"Read:   {read - start:.3f}s")
print(f"Create: {created - read:.3f}s")
print(f"Write:  {done - created:.3f}s")
print(f"Total:  {done - start:.3f}s")

Read:   2.320s
Create: 43.147s
Write:  7.384s
Total:  52.851s

PDAL

Code

import time

from pdal import Filter, Reader

pipeline = Reader.las(
    filename="USGS_LPC_CO_SoPlatteRiver_Lot2a_2013_13TDE489446_LAS_2015.copc.laz"
) | Filter.range(limits="Classification[2:2]")

start = time.time()
pipeline.execute()

read = time.time()
pipeline = Filter.delaunay().pipeline(pipeline.arrays[0])
pipeline.execute()

created = time.time()
mesh = pipeline.get_meshio(0)
mesh.write("dt-pdal.ply")

done = time.time()

print(f"Read:   {read - start:.3f}s")
print(f"Create: {created - read:.3f}s")
print(f"Write:  {done - created:.3f}s")
print(f"Total:  {done - start:.3f}s")

Read:   3.916s
Create: 3.943s
Write:  1.128s
Total:  8.987s

QC

Comparing statistics

$ pdal info dt-startinpy.ply | jq .stats.statistic 
[
  {
    "average": 489726.2401,
    "count": 2872483,
    "maximum": 490500,
    "minimum": 489000,
    "name": "X",
    "position": 0,
    "stddev": 432.9948745,
    "variance": 187484.5614
  },
  {
    "average": 4446745.876,
    "count": 2872483,
    "maximum": 4447500,
    "minimum": 4446000,
    "name": "Y",
    "position": 1,
    "stddev": 439.3933042,
    "variance": 193066.4758
  },
  {
    "average": 1533.628089,
    "count": 2872483,
    "maximum": 1554.64,
    "minimum": 1515.02,
    "name": "Z",
    "position": 2,
    "stddev": 7.959375129,
    "variance": 63.35165245
  }
]
$ pdal info dt-pdal.ply | jq .stats.statistic
[
  {
    "average": 489726.2391,
    "count": 2872606,
    "maximum": 490500,
    "minimum": 489000,
    "name": "X",
    "position": 0,
    "stddev": 432.996162,
    "variance": 187485.6763
  },
  {
    "average": 4446745.877,
    "count": 2872606,
    "maximum": 4447500,
    "minimum": 4446000,
    "name": "Y",
    "position": 1,
    "stddev": 439.392411,
    "variance": 193065.6908
  },
  {
    "average": 1533.628108,
    "count": 2872606,
    "maximum": 1554.64,
    "minimum": 1515.02,
    "name": "Z",
    "position": 2,
    "stddev": 7.959359551,
    "variance": 63.35140446
  }
]

Images from QGIS

The input point cloud:

The output mesh (both startinpy and PDAL produced the same-enough output):

[hobu]

Filter.range(limits="Classification[2:2]")

Filters.expression(expression="Classification == 2") <-- much better modern PDAL way to the same thing.

[hugoledoux]

added : b33d1cb