Refine new point cloud gridding function with comments (#558)

geoutils/pointcloud.py

@@ -12,7 +12,7 @@
 
 def _grid_pointcloud(
     pc: gpd.GeoDataFrame,
-    grid_coords: tuple[NDArrayNum, NDArrayNum],
+    grid_coords: tuple[NDArrayNum, NDArrayNum] = None,
     data_column_name: str = "b1",
     resampling: Literal["nearest", "linear", "cubic"] = "linear",
     dist_nodata_pixel: float = 1.0,
@@ -25,13 +25,27 @@ def _grid_pointcloud(
     matter the distance).
 
     :param pc: Point cloud.
-    :param grid_coords: Grid coordinates for X and Y.
+    :param grid_coords: Regular raster grid coordinates in X and Y (i.e. equally spaced, independently for each axis).
     :param data_column_name: Name of data column for point cloud (only if passed as a geodataframe).
     :param resampling: Resampling method within delauney triangles (defaults to linear).
     :param dist_nodata_pixel: Distance from the point cloud after which grid cells are filled by nodata values,
         expressed in number of pixels.
     """
 
+    # Input checks
+    if (
+        not isinstance(grid_coords, tuple)
+        or not (isinstance(grid_coords[0], np.ndarray) and grid_coords[0].ndim == 1)
+        or not (isinstance(grid_coords[1], np.ndarray) and grid_coords[1].ndim == 1)
+    ):
+        raise TypeError("Input grid coordinates must be 1D arrays.")
+
+    diff_x = np.diff(grid_coords[0])
+    diff_y = np.diff(grid_coords[1])
+
+    if not all(diff_x == diff_x[0]) and all(diff_y == diff_y[0]):
+        raise ValueError("Grid coordinates must be regular (equally spaced, independently along X and Y).")
+
     # 1/ Interpolate irregular point cloud on a regular grid
 
     # Get meshgrid coordinates

tests/test_pointcloud.py

@@ -2,6 +2,7 @@
 
 import geopandas as gpd
 import numpy as np
+import pytest
 import rasterio as rio
 from shapely import geometry
 
@@ -10,7 +11,7 @@
 
 
 class TestPointCloud:
-    def test_grid_pc__chull(self) -> None:
+    def test_grid_pc(self) -> None:
         """Test point cloud gridding."""
 
         # 1/ Check gridding interpolation falls back exactly on original raster
@@ -46,7 +47,7 @@ def test_grid_pc__chull(self) -> None:
         poly = geometry.MultiPoint([[p.x, p.y] for p in pc.geometry])
         chull = poly.convex_hull
 
-        # We compute the index of grid cells interesting the convex hull
+        # We compute the index of grid cells intersecting the convex hull
         ind_inters_convhull = rst_pc.intersects(chull)
 
         # We get corresponding 1D indexes for gridded output
@@ -105,3 +106,10 @@ def test_grid_pc__chull(self) -> None:
         ifarchull, jfarchull = list(zip(*far_in_chull))
 
         assert all(~np.isfinite(gridded_pc[ifarchull, jfarchull]))
+
+        # 3/ Errors
+        with pytest.raises(TypeError, match="Input grid coordinates must be 1D arrays.*"):
+            Raster.from_pointcloud_regular(pc, grid_coords=(1, "lol"))  # type: ignore
+        with pytest.raises(ValueError, match="Grid coordinates must be regular*"):
+            grid_coords[0][0] += 1
+            Raster.from_pointcloud_regular(pc, grid_coords=grid_coords)  # type: ignore