Add PointCloud subclass of Vector

dev-environment.yml

@@ -13,6 +13,7 @@ dependencies:
   - tqdm
   - xarray
   - rioxarray=0.*
+  - python-pdal
 
   # Development-specific, to mirror manually in setup.cfg [options.extras_require].
   - pip

environment.yml

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

geoutils/vector.py

@@ -26,6 +26,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
+import pdal
 import rasterio as rio
 import rasterio.errors
 import shapely
@@ -54,7 +55,7 @@
 
 class Vector:
     """
-    The georeferenced vector
+    The georeferenced vector.
 
      Main attributes:
         ds: :class:`geopandas.GeoDataFrame`
@@ -1665,6 +1666,129 @@ def buffer_without_overlap(self, buffer_size: int | float, metric: bool = True,
         return Vector(merged_voronoi)
 
 
+######################
+# Point cloud subclass
+######################
+
+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(Vector):
+    """
+    The georeferenced point cloud.
+
+    A point cloud is a vector of 2D point geometries associated to values from a specific data column.
+
+     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 Raster 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
+
+    @classmethod
+    def from_array(cls, array: np.ndarray, crs: CRS, data_column: str | None = None) -> 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={"z": 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)."""
+
+        array_in = np.array(zip(*tuples))
+
+
+    def to_array(self):
+        """Convert point cloud to a 3 x N array of X coordinates, Y coordinates and Z values."""
+
+
+        return
+
+    def to_tuples(self):
+        """Convert point cloud to a N-size list of tuples (X coordinate, Y coordinate, Z value)."""
+
+        return
+
+    def interp_raster(self, raster: gu.Raster):
+        """Interpolate raster at point cloud coordinates."""
+
+        return
+
 # -----------------------------------------
 # Additional stand-alone utility functions
 # -----------------------------------------