✨ LaspyReader for reading LAS/LAZ point cloud files

docs/_config.yml

@@ -54,6 +54,9 @@ sphinx:
       geopandas:
         - 'https://geopandas.org/en/latest/'
         - null
+      laspy:
+        - 'https://laspy.readthedocs.io/en/latest/'
+        - null
       mmdetection:
         - 'https://mmdetection.readthedocs.io/zh_CN/latest/'
         - null

docs/api.md

@@ -27,6 +27,15 @@
     :show-inheritance:
 ```
 
+### Laspy
+
+```{eval-rst}
+.. automodule:: zen3geo.datapipes.laspy
+.. autoclass:: zen3geo.datapipes.LaspyReader
+.. autoclass:: zen3geo.datapipes.laspy.LaspyReaderIterDataPipe
+    :show-inheritance:
+```
+
 ### Pyogrio
 
 ```{eval-rst}

docs/index.md

@@ -7,14 +7,15 @@ Get what you need, not more, not less:
 | Command                        |  Dependencies |
 |:-------------------------------|---------------|
 | `pip install zen3geo`          | rioxarray, torchdata |
+| `pip install zen3geo[lidar]`   | rioxarray, torchdata, laspy[lazrs] |
 | `pip install zen3geo[raster]`  | rioxarray, torchdata, xbatcher, zarr |
 | `pip install zen3geo[spatial]` | rioxarray, torchdata, datashader, spatialpandas |
 | `pip install zen3geo[stac]`    | rioxarray, torchdata, pystac, pystac-client, stackstac, xpystac |
 | `pip install zen3geo[vector]`  | rioxarray, torchdata, pyogrio[geopandas] |
 
 Retrieve more ['extras'](https://github.com/weiji14/zen3geo/blob/main/pyproject.toml) using
 
-    pip install zen3geo[raster,spatial,stac,vector]
+    pip install zen3geo[lidar,raster,spatial,stac,vector]
 
 To install the development version from [TestPyPI](https://test.pypi.org/project/zen3geo), do:
 

pyproject.toml

@@ -35,6 +35,7 @@ rioxarray = ">=0.10.0"
 torchdata = ">=0.4.0"
 # Optional
 datashader = {version = ">=0.14.0", optional = true}
+laspy = {version = ">=2.5.0", extras = ["lazrs"], optional = true}
 pyogrio = {version = ">=0.4.0", extras = ["geopandas"], optional = true}
 pystac = {version=">=1.4.0", optional=true}
 pystac-client = {version = ">=0.4.0", optional = true}
@@ -76,6 +77,9 @@ docs = [
     "xpystac",
     "zarr"
 ]
+lidar = [
+    "laspy",
+]
 raster = [
     "xbatcher",
     "zarr"

zen3geo/datapipes/__init__.py

@@ -9,6 +9,7 @@
 from zen3geo.datapipes.geopandas import (
     GeoPandasRectangleClipperIterDataPipe as GeoPandasRectangleClipper,
 )
+from zen3geo.datapipes.laspy import LaspyReaderIterDataPipe as LaspyReader
 from zen3geo.datapipes.pyogrio import PyogrioReaderIterDataPipe as PyogrioReader
 from zen3geo.datapipes.pystac import PySTACItemReaderIterDataPipe as PySTACItemReader
 from zen3geo.datapipes.pystac_client import (

zen3geo/datapipes/laspy.py

@@ -0,0 +1,89 @@
+"""
+DataPipes for :doc:`laspy <laspy:index>`.
+"""
+from typing import Any, Dict, Iterator, Optional
+
+try:
+    import laspy
+except ImportError:
+    laspy = None
+from torchdata.datapipes import functional_datapipe
+from torchdata.datapipes.iter import IterDataPipe
+from torchdata.datapipes.utils import StreamWrapper
+
+
+@functional_datapipe("read_from_laspy")
+class LaspyReaderIterDataPipe(IterDataPipe[StreamWrapper]):
+    """
+    Takes LAS/LAZ files from local disk or an :py:class:`io.BytesIO` stream (as long as
+    they can be read by laspy) and yields :py:class:`laspy.lasdata.LasData` objects
+    (functional name: ``read_from_laspy``).
+
+    Parameters
+    ----------
+    source_datapipe : IterDataPipe[str, io.BytesIO]
+        A DataPipe that contains filepaths or an :py:class:`io.BytesIO` stream to point
+        cloud data such as LAS, LAZ, COPC, etc.
+
+    kwargs : Optional
+        Extra keyword arguments to pass to :py:func:`laspy.read`.
+
+    Yields
+    ------
+    stream_obj : laspy.lasdata.LasData
+        A :py:class:`laspy.lasdata.LasData` object containing the point cloud data.
+
+    Raises
+    ------
+    ModuleNotFoundError
+        If ``laspy`` is not installed. See
+        :doc:`install instructions for laspy <laspy:installation>`, (e.g. via
+        ``pip install laspy[lazrs]``) before using this class.
+
+    Example
+    -------
+    >>> import pytest
+    >>> laspy = pytest.importorskip("laspy")
+    ...
+    >>> from torchdata.datapipes.iter import IterableWrapper
+    >>> from zen3geo.datapipes import LaspyReader
+    ...
+    >>> # Read in LAZ data using DataPipe
+    >>> file_url: str = "https://opentopography.s3.sdsc.edu/pc-bulk/NZ19_Wellington/CL2_BQ31_2019_1000_2138.laz"
+    >>> dp = IterableWrapper(iterable=[file_url])
+    >>> _, dp_stream = dp.read_from_http().unzip(sequence_length=2)
+    >>> dp_laspy = dp_stream.read_from_laspy()
+    ...
+    >>> # Loop or iterate over the DataPipe stream
+    >>> it = iter(dp_laspy)
+    >>> lasdata = next(it)
+    >>> lasdata.header
+    <LasHeader(1.4, <PointFormat(6, 0 bytes of extra dims)>)>
+    >>> lasdata.xyz
+    array([[ 1.74977156e+06,  5.42749877e+06, -7.24000000e-01],
+           [ 1.74977152e+06,  5.42749846e+06, -7.08000000e-01],
+           [ 1.74977148e+06,  5.42749815e+06, -7.00000000e-01],
+           ...,
+           [ 1.74976026e+06,  5.42756798e+06, -4.42000000e-01],
+           [ 1.74976029e+06,  5.42756829e+06, -4.17000000e-01],
+           [ 1.74976032e+06,  5.42756862e+06, -4.04000000e-01]])
+    """
+
+    def __init__(
+        self, source_datapipe: IterDataPipe[str], **kwargs: Optional[Dict[str, Any]]
+    ) -> None:
+        if laspy is None:
+            raise ModuleNotFoundError(
+                "Package `laspy` is required to be installed to use this datapipe. "
+                "Please use `pip install laspy` or "
+                "`conda install -c conda-forge laspy` to install the package"
+            )
+        self.source_datapipe: IterDataPipe[str] = source_datapipe
+        self.kwargs = kwargs
+
+    def __iter__(self) -> Iterator[StreamWrapper]:
+        for lazstream in self.source_datapipe:
+            yield StreamWrapper(laspy.read(source=lazstream, **self.kwargs))
+
+    def __len__(self) -> int:
+        return len(self.source_datapipe)

zen3geo/tests/test_datapipes_laspy.py

@@ -0,0 +1,119 @@
+"""
+Tests for laspy datapipes.
+"""
+import tempfile
+import urllib
+
+import numpy as np
+import numpy.testing as npt
+import pytest
+from torchdata.datapipes.iter import IterableWrapper
+
+from zen3geo.datapipes import LaspyReader
+
+laspy = pytest.importorskip("laspy")
+
+
+# %%
+def test_laspy_reader_las_local():
+    """
+    Ensure that LaspyReader works to read in a LAS file (on disk) and outputs a
+    laspy.lasdata.LasData object.
+    """
+    with tempfile.NamedTemporaryFile(suffix=".las") as tmpfile:
+        urllib.request.urlretrieve(
+            url="https://github.com/laz-rs/laz-rs/raw/0.8.3/tests/data/point-time-color.las",
+            filename=tmpfile.name,
+        )
+        dp = IterableWrapper(iterable=[tmpfile.name])
+
+        # Using class constructors
+        dp_laspy = LaspyReader(source_datapipe=dp)
+        # Using functional form (recommended)
+        dp_laspy = dp.read_from_laspy()
+
+        assert len(dp_laspy) == 1
+        it = iter(dp_laspy)
+        lasdata = next(it)
+
+    assert lasdata.header.version == laspy.header.Version(major=1, minor=2)
+    assert lasdata.header.point_format == laspy.point.PointFormat(point_format_id=3)
+    assert lasdata.points.array.shape == (1065,)
+    assert lasdata.xyz.shape == (1065, 3)
+    npt.assert_allclose(
+        actual=lasdata.xyz.mean(axis=0),
+        desired=[494494.6635117371, 4878134.831230047, 132.31299530516432],
+    )
+    npt.assert_allclose(actual=np.unique(lasdata.classification.array), desired=[1, 2])
+
+
+def test_laspy_reader_laz_http():
+    """
+    Ensure that LaspyReader works to read in a LAZ file (from a HTTP byte stream)
+    and outputs a laspy.lasdata.LasData object.
+    """
+    file_url: str = "https://github.com/laz-rs/laz-rs/raw/0.8.3/tests/data/point-version-1-point-wise.laz"
+    dp = IterableWrapper(iterable=[file_url])
+    _, dp_stream = (
+        dp.read_from_http()
+        .read_from_stream()
+        .set_length(length=1)
+        .unzip(sequence_length=2)
+    )
+
+    # Using class constructors
+    dp_laspy = LaspyReader(source_datapipe=dp_stream)
+    # Using functional form (recommended)
+    dp_laspy = dp_stream.read_from_laspy()
+
+    assert len(dp_laspy) == 1
+    it = iter(dp_laspy)
+    lasdata = next(it)
+
+    assert lasdata.header.version == laspy.header.Version(major=1, minor=0)
+    assert lasdata.header.point_format == laspy.point.PointFormat(point_format_id=0)
+    assert lasdata.points.array.shape == (11781,)
+    assert lasdata.xyz.shape == (11781, 3)
+    npt.assert_allclose(
+        actual=lasdata.xyz.mean(axis=0),
+        desired=[2483799.026934895, 366405.56612511666, 1511.9428214922332],
+    )
+    npt.assert_allclose(
+        actual=np.unique(lasdata.classification), desired=[1, 2, 8, 9, 12, 15]
+    )
+
+
+def test_laspy_reader_copc_http():
+    """
+    Ensure that LaspyReader works to read in a COPC file (from a HTTP byte stream) and
+    outputs a laspy.lasdata.LasData object.
+    """
+    file_url: str = (
+        "https://github.com/laspy/laspy/raw/2.5.3/tests/data/simple_with_page.copc.laz"
+    )
+    dp = IterableWrapper(iterable=[file_url])
+    _, dp_stream = (
+        dp.read_from_http()
+        .read_from_stream()
+        .set_length(length=1)
+        .unzip(sequence_length=2)
+    )
+
+    # Using class constructors
+    dp_laspy = LaspyReader(source_datapipe=dp_stream)
+    # Using functional form (recommended)
+    dp_laspy = dp_stream.read_from_laspy()
+
+    assert len(dp_laspy) == 1
+    it = iter(dp_laspy)
+    lasdata = next(it)
+
+    assert lasdata.header.version == laspy.header.Version(major=1, minor=4)
+    assert lasdata.header.point_format == laspy.point.PointFormat(point_format_id=7)
+    assert lasdata.points.array.shape == (1065,)
+    assert lasdata.xyz.shape == (1065, 3)
+    npt.assert_allclose(
+        actual=lasdata.xyz.mean(axis=0),
+        desired=[637296.7351830985, 851249.5384882629, 434.0978403755869],
+    )
+    npt.assert_allclose(actual=np.unique(lasdata.classification), desired=[1, 2])