Incremental commit on PointCloud class

dev-environment.yml

@@ -14,12 +14,12 @@ dependencies:
   - xarray
   - dask
   - rioxarray=0.*
-  - python-pdal
 
   # Development-specific, to mirror manually in setup.cfg [options.extras_require].
   - pip
 
   # Optional dependencies
+  - laspy
   - scikit-image
 
   # Test dependencies
@@ -50,4 +50,4 @@ dependencies:
   - typing-extensions
 
   - pip:
-    - -e ./
+      - -e ./

environment.yml

@@ -14,4 +14,3 @@ dependencies:
   - xarray
   - dask
   - rioxarray=0.*
-  - python-pdal

geoutils/__init__.py

@@ -2,10 +2,11 @@
 GeoUtils is a Python package for the analysis of geospatial data.
 """
 
-from geoutils import examples, pointcloud, projtools, raster, vector  # noqa
+from geoutils import examples  # noqa
 from geoutils._config import config  # noqa
 from geoutils.raster import Mask, Raster  # noqa
 from geoutils.vector import Vector  # noqa
+from geoutils.pointcloud import PointCloud  # noqa
 
 try:
     from geoutils._version import __version__ as __version__  # noqa

geoutils/interface/gridding.py

@@ -3,8 +3,10 @@
 import warnings
 from typing import Literal
 
+import affine
 import geopandas as gpd
 import numpy as np
+import rasterio as rio
 from scipy.interpolate import griddata
 
 from geoutils._typing import NDArrayNum
@@ -16,7 +18,7 @@ def _grid_pointcloud(
     data_column_name: str = "b1",
     resampling: Literal["nearest", "linear", "cubic"] = "linear",
     dist_nodata_pixel: float = 1.0,
-) -> NDArrayNum:
+) -> tuple[NDArrayNum, affine.Affine]:
     """
     Grid point cloud (possibly irregular coordinates) to raster (regular grid) using delaunay triangles interpolation.
 
@@ -88,4 +90,7 @@ def _grid_pointcloud(
     # Flip Y axis of grid
     aligned_dem = np.flip(aligned_dem, axis=0)
 
-    return aligned_dem
+    # 3/ Derive output transform from input grid
+    transform_from_coords = rio.transform.from_origin(grid_coords[0][0], grid_coords[0][1], res_x, res_y)
+
+    return aligned_dem, transform_from_coords

geoutils/pointcloud/__init__.py

@@ -1 +1 @@
-from geoutils.pointcloud.pointcloud import *  # noqa
+from geoutils.pointcloud.pointcloud import PointCloud  # noqa

geoutils/pointcloud/pointcloud.py

@@ -1 +1,179 @@
-"""Module for future PointCloud class."""
+"""Module for PointCloud class."""
+from __future__ import annotations
+
+import warnings
+from typing import Any, Iterable, Literal
+import pathlib
+
+from pyproj import CRS
+import numpy as np
+import pandas as pd
+import geopandas as gpd
+from shapely.geometry.base import BaseGeometry
+
+import geoutils as gu
+from geoutils.interface.gridding import _grid_pointcloud
+from geoutils._typing import Number
+
+try:
+    import laspy
+    has_laspy = True
+except ImportError:
+    has_laspy = False
+
+
+def _read_pdal(filename: str, **kwargs: Any) -> gpd.GeoDataFrame:
+    """Read a point cloud dataset with PDAL."""
+
+    # Basic json command to read an entire file
+    json_string = f"""
+    [
+        "{filename}"
+    ]
+    """
+
+    # Run and extract array as dataframe
+    pipeline = pdal.Pipeline(json_string)
+    pipeline.execute()
+    df = pd.DataFrame(pipeline.arrays[0])
+
+    # Build geodataframe from 2D points and data table
+    crs = CRS.from_wkt(pipeline.srswkt2)
+    gdf = gpd.GeoDataFrame(geometry=gpd.points_from_xy(x=df["X"], y=df["Y"], crs=crs), data=df.iloc[:, 2:])
+
+    return gdf
+
+def _write_pdal(filename: str, **kwargs):
+    """Write a point cloud dataset with PDAL."""
+
+    return
+
+
+class PointCloud(gu.Vector):
+    """
+    The georeferenced point cloud.
+
+    A point cloud is a vector of 2D point geometries associated to values from a main data column, with access
+    to values from auxiliary data columns.
+
+     Main attributes:
+        ds: :class:`geopandas.GeoDataFrame`
+            Geodataframe of the point cloud.
+        data_column: str
+            Name of point cloud data column.
+        crs: :class:`pyproj.crs.CRS`
+            Coordinate reference system of the point cloud.
+        bounds: :class:`rio.coords.BoundingBox`
+            Coordinate bounds of the point cloud.
+
+
+    All other attributes are derivatives of those attributes, or read from the file on disk.
+    See the API for more details.
+    """
+
+    _data_column: str | None
+
+    def __init__(self,
+                 filename_or_dataset: str | pathlib.Path | gpd.GeoDataFrame | gpd.GeoSeries | BaseGeometry,
+                 data_column: str | None = None,):
+        """
+        Instantiate a point cloud from either a data column name and a vector (filename, GeoPandas dataframe or series,
+        or a Shapely geometry), or only with a point cloud file type.
+
+        :param filename_or_dataset: Path to file, or GeoPandas dataframe or series, or Shapely geometry.
+        :param data_column: Name of data column defining the point cloud.
+        """
+
+        # If PointCloud is passed, simply point back to PointCloud
+        if isinstance(filename_or_dataset, PointCloud):
+            for key in filename_or_dataset.__dict__:
+                setattr(self, key, filename_or_dataset.__dict__[key])
+            return
+        # Else rely on parent Vector class options (including raised errors)
+        else:
+            super().__init__(filename_or_dataset)
+
+        if data_column not in self.ds.columns():
+            raise ValueError(f"Data column {data_column} not found in vector file, available columns "
+                             f"are: {', '.join(self.ds.columns)}.")
+        self._data_column = data_column
+
+    @property
+    def data(self):
+        return self.ds[self.data_column].values
+
+    @property
+    def data_column(self) -> str:
+        """Name of main data column of the point cloud."""
+        return self._data_column
+
+    @property
+    def auxiliary_columns(self) -> list[str]:
+        """Name of auxiliary data columns of the point cloud."""
+        return self._auxiliary_columns
+
+    @classmethod
+    def from_array(cls, array: np.ndarray, crs: CRS, data_column: str | None = "z") -> PointCloud:
+        """Create point cloud from a 3 x N or N x 3 array of X coordinate, Y coordinates and Z values."""
+
+        # Check shape
+        if array.shape[0] != 3 and array.shape[1] != 3:
+            raise ValueError("Array must be of shape 3xN or Nx3.")
+
+        # Make the first axis the one with size 3
+        if array.shape[0] != 3:
+            array_in = array.T
+        else:
+            array_in = array
+
+        # Build geodataframe
+        gdf = gpd.GeoDataFrame(geometry=gpd.points_from_xy(x=array_in[0, :], y=array_in[1, :], crs=crs),
+                               data={data_column: array[2, :]})
+
+        return cls(filename_or_dataset=gdf, data_column=data_column)
+
+
+    @classmethod
+    def from_tuples(self, tuples: Iterable[tuple[Number, Number, Number]], crs: CRS, data_column: str | None = None):
+        """Create point cloud from a N-size list of tuples (X coordinate, Y coordinate, Z value)."""
+
+        self.from_array(np.array(zip(*tuples)), crs=crs, data_column=data_column)
+
+
+    def to_array(self):
+        """Convert point cloud to a 3 x N array of X coordinates, Y coordinates and Z values."""
+
+        return np.stack((self.geometry.x.values, self.geometry.y.values, self.ds[self.data_column].values), axis=0)
+
+    def to_tuples(self):
+        """Convert point cloud to a N-size list of tuples (X coordinate, Y coordinate, Z value)."""
+
+        return self.geometry.x.values, self.geometry.y.values, self.ds[self.data_column].values
+
+    def grid(self,
+             ref: gu.Raster | None,
+             grid_coords: tuple[np.ndarray, np.ndarray] | None,
+             resampling: Literal["nearest", "linear", "cubic"],
+             dist_nodata_pixel: float = 1.) -> gu.Raster:
+        """Grid point cloud into a raster."""
+
+        if isinstance(ref, gu.Raster):
+            if grid_coords is None:
+                warnings.warn("Both reference raster and grid coordinates were passed for gridding, "
+                              "using only the reference raster.")
+            grid_coords = ref.coords(grid=False)
+        else:
+            grid_coords = grid_coords
+
+        array, transform = _grid_pointcloud(self.ds, grid_coords=grid_coords, data_column_name=self.data_column,
+                                 resampling=resampling, dist_nodata_pixel=dist_nodata_pixel)
+
+        return gu.Raster.from_array(data=array, transform=transform, crs=self.crs, nodata=None)
+
+    # def subsample(self) -> PointCloud | NDArrayf:
+    #
+    # @classmethod
+    # def from_raster(cls, raster: gu.Raster) -> PointCloud:
+    #     """Create a point cloud from a raster. Equivalent with Raster.to_pointcloud."""
+    #
+    #    pc = _raster_to_pointcloud(source_raster=raster)

tests/test_interface/test_gridding.py

@@ -1,4 +1,4 @@
-"""Test module for point cloud functionalities."""
+"""Test module for gridding functionalities."""
 
 import geopandas as gpd
 import numpy as np
@@ -30,11 +30,13 @@ def test_grid_pc(self) -> None:
         grid_coords = rst.coords(grid=False)
 
         # Grid the point cloud
-        gridded_pc = _grid_pointcloud(pc, grid_coords=grid_coords)
+        gridded_pc, output_transform = _grid_pointcloud(pc, grid_coords=grid_coords)
 
         # Compare back to raster, all should be very close (but not exact, some info is lost due to interpolations)
         valids = np.isfinite(gridded_pc)
         assert np.allclose(gridded_pc[valids], rst.data.data[valids], rtol=10e-5)
+        # And the transform exactly the same
+        assert output_transform == transform
 
         # 2/ Check the propagation of nodata values
 
@@ -86,7 +88,7 @@ def test_grid_pc(self) -> None:
         assert all(~np.isfinite(gridded_pc[ifarchull, jfarchull]))
 
         # Check for a different distance value
-        gridded_pc = _grid_pointcloud(pc, grid_coords=grid_coords, dist_nodata_pixel=0.5)
+        gridded_pc, output_transform = _grid_pointcloud(pc, grid_coords=grid_coords, dist_nodata_pixel=0.5)
         ind_close = np.array(list_min_dist) <= 0.5
 
         # We get the indexes for these coordinates

tests/test_pointcloud/test_pointcloud.py

@@ -1 +1,17 @@
-"""Test for future PointCloud class."""
+"""Test for PointCloud class."""
+
+class TestPointCloud:
+
+    def test_init(self):
+
+    def test_data_column_name(self):
+
+    def test_auxiliary_column_name(self):
+
+    def test_from_array(self):
+
+    def test_from_tuples(self):
+
+    def test_to_array(self):
+
+    def test_to_tuples(self):