Small fixes and optimizations for to_pointcloud and interp_points (…

…#549)

doc/source/raster_class.md

@@ -317,7 +317,7 @@ rast_reproj.value_at_coords(x=0.5, y=0.5, window=3, reducer_function=np.ma.media
 
 ```{code-cell} ipython3
 # Interpolate coordinate value with quintic algorithm
-rast_reproj.interp_points([(0.5, 0.5)], method="quintic")
+rast_reproj.interp_points((0.5, 0.5), method="quintic")
 ```
 
 ```{note}

examples/analysis/point_extraction/interpolation.py

@@ -28,7 +28,7 @@
 x_coords = rng.uniform(rast.bounds.left + 50, rast.bounds.right - 50, 50)
 y_coords = rng.uniform(rast.bounds.bottom + 50, rast.bounds.top - 50, 50)
 
-vals = rast.interp_points(points=list(zip(x_coords, y_coords)))
+vals = rast.interp_points(points=(x_coords, y_coords))
 
 # %%
 # Replace by Vector function once done
@@ -50,7 +50,7 @@
 # %%
 # We can interpolate again by shifting according to our interpretation, and changing the resampling algorithm (default to "linear").
 
-vals_shifted = rast.interp_points(points=list(zip(x_coords, y_coords)), shift_area_or_point=True, method="quintic")
+vals_shifted = rast.interp_points(points=(x_coords, y_coords), shift_area_or_point=True, method="quintic")
 np.nanmean(vals - vals_shifted)
 
 # %%

geoutils/raster/raster.py

@@ -3657,7 +3657,7 @@ def outside_image(self, xi: ArrayLike, yj: ArrayLike, index: bool = True) -> boo
 
     def interp_points(
         self,
-        points: tuple[list[float], list[float]],
+        points: tuple[Number, Number] | tuple[NDArrayNum, NDArrayNum],
         method: Literal["nearest", "linear", "cubic", "quintic"] = "linear",
         band: int = 1,
         input_latlon: bool = False,
@@ -3675,8 +3675,8 @@ def interp_points(
          to ensure that the interpolation of points is done at the right location. See parameter description
          of shift_area_or_point for more details.
 
-        :param points: Point(s) at which to interpolate raster value. If points fall outside of image, value
-            returned is nan. Shape should be (N,2).
+        :param points: Point(s) at which to interpolate raster value (tuple of X/Y array-likes). If points fall
+        outside of image, value returned is nan.
         :param method: Interpolation method, one of 'nearest', 'linear', 'cubic', or 'quintic'. For more information,
             see scipy.ndimage.map_coordinates and scipy.interpolate.interpn. Default is linear.
         :param band: Band to use (from 1 to self.count).
@@ -3690,7 +3690,7 @@ def interp_points(
         """
 
         # Get coordinates
-        x, y = list(zip(*points))
+        x, y = points
 
         # If those are in latlon, convert to Raster CRS
         if input_latlon:
@@ -3953,6 +3953,10 @@ def to_pointcloud(
             all_bands = [data_band]
             all_column_names = [data_column_name]
 
+        # If subsample is the entire array, load it to optimize speed
+        if subsample == 1 and not self.is_loaded:
+            self.load(bands=all_bands)
+
         # Band indexes in the array are band number minus one
         all_indexes = [b - 1 for b in all_bands]
 
@@ -4008,8 +4012,9 @@ def to_pointcloud(
         if np.ma.isMaskedArray(pixel_data):
             pixel_data = pixel_data.data
 
-        # If nodata values were not skipped, convert them to NaNs
+        # If nodata values were not skipped, convert them to NaNs and change data type
         if skip_nodata is False:
+            pixel_data = pixel_data.astype("float32")
             pixel_data[pixel_data == self.nodata] = np.nan
 
         # Now we force the coordinates we define for the point cloud, according to pixel interpretation
@@ -4048,7 +4053,7 @@ def from_pointcloud_regular(
         Create a raster from a point cloud with coordinates on a regular grid.
 
         To inform on what grid to create the raster, either pass a tuple of X/Y grid coordinates, or the expected
-        transform and shape. All point cloud coordinates must fall exactly at one the coordinates of this grid.
+        transform and shape. All point cloud coordinates must fall exactly at one of the coordinates of this grid.
 
         :param pointcloud: Point cloud.
         :param grid_coords: Regular coordinate vectors for the raster, from which the geotransform and shape are

tests/test_raster/test_raster.py

@@ -2070,7 +2070,6 @@ def test_xy2ij(self) -> None:
         rng = np.random.default_rng(42)
         xrand = rng.integers(low=0, high=r.width, size=(10,)) * list(r.transform)[0] + xmin
         yrand = ymax + rng.integers(low=0, high=r.height, size=(10,)) * list(r.transform)[4]
-        pts = list(zip(xrand, yrand))
 
         # Get decimal indexes based on "Point", should refer to the corner still (shift False by default)
         i, j = r.xy2ij(xrand, yrand, shift_area_or_point=False)
@@ -2102,15 +2101,14 @@ def test_xy2ij(self) -> None:
             list_z_ind.append(z_ind)
 
         # First order interpolation
-        rpts = r.interp_points(pts, method="linear")
+        rpts = r.interp_points((xrand, yrand), method="linear")
         # The values interpolated should be equal
         assert np.array_equal(np.array(list_z_ind, dtype=np.float32), rpts, equal_nan=True)
 
         # Test there is no failure with random coordinates (edge effects, etc)
         xrand = rng.uniform(low=xmin, high=xmax, size=(1000,))
         yrand = rng.uniform(low=ymin, high=ymax, size=(1000,))
-        pts = list(zip(xrand, yrand))
-        r.interp_points(pts)
+        r.interp_points((xrand, yrand))
 
         # Second, test after a crop: the Raster now has an Area interpretation, those should fall right on the integer
         # pixel indexes
@@ -2124,7 +2122,6 @@ def test_xy2ij(self) -> None:
         # give back the same values that fall right on the coordinates
         xrand = rng.integers(low=0, high=r.width, size=(10,)) * list(r.transform)[0] + xmin
         yrand = ymax + rng.integers(low=0, high=r.height, size=(10,)) * list(r.transform)[4]
-        pts = list(zip(xrand, yrand))
         # By default, i and j are returned as integers
         i, j = r.xy2ij(xrand, yrand, op=np.float32)
         list_z_ind = []
@@ -2134,21 +2131,21 @@ def test_xy2ij(self) -> None:
             z_ind = img[int(i[k]), int(j[k])]
             list_z_ind.append(z_ind)
 
-        rpts = r.interp_points(pts, method="linear")
+        rpts = r.interp_points((xrand, yrand), method="linear")
 
         assert np.array_equal(np.array(list_z_ind, dtype=np.float32), rpts, equal_nan=True)
 
         # Test for an invidiual point (shape can be tricky at 1 dimension)
         x = 493120.0
         y = 3101000.0
         i, j = r.xy2ij(x, y)
-        val = r.interp_points([(x, y)], method="linear")[0]
+        val = r.interp_points((x, y), method="linear")[0]
         val_img = img[int(i[0]), int(j[0])]
         assert val_img == val
 
         # Finally, check that interp convert to latlon
         lat, lon = gu.projtools.reproject_to_latlon([x, y], in_crs=r.crs)
-        val_latlon = r.interp_points([(lat, lon)], method="linear", input_latlon=True)[0]
+        val_latlon = r.interp_points((lat, lon), method="linear", input_latlon=True)[0]
         assert val == pytest.approx(val_latlon, abs=0.0001)
 
     @pytest.mark.parametrize("tag_aop", [None, "Area", "Point"])  # type: ignore
@@ -2171,19 +2168,17 @@ def test_interp_points__synthetic(self, tag_aop: str | None, shift_aop: bool) ->
         # The actual X/Y coords will be offset by one because Y axis is inverted and pixel coords is upper-left corner
         points_x, points_y = raster.ij2xy(i=index_x, j=index_y, shift_area_or_point=shift_aop)
 
-        points = np.array((points_x, points_y)).T
-
         # The following 4 methods should yield the same result because:
         # Nearest = Linear interpolation at the location of a data point
         # Regular grid = Equal grid interpolation at the location of a data point
 
-        raster_points = raster.interp_points(points, method="nearest", shift_area_or_point=shift_aop)
-        raster_points_lin = raster.interp_points(points, method="linear", shift_area_or_point=shift_aop)
+        raster_points = raster.interp_points((points_x, points_y), method="nearest", shift_area_or_point=shift_aop)
+        raster_points_lin = raster.interp_points((points_x, points_y), method="linear", shift_area_or_point=shift_aop)
         raster_points_interpn = raster.interp_points(
-            points, method="nearest", force_scipy_function="interpn", shift_area_or_point=shift_aop
+            (points_x, points_y), method="nearest", force_scipy_function="interpn", shift_area_or_point=shift_aop
         )
         raster_points_interpn_lin = raster.interp_points(
-            points, method="linear", force_scipy_function="interpn", shift_area_or_point=shift_aop
+            (points_x, points_y), method="linear", force_scipy_function="interpn", shift_area_or_point=shift_aop
         )
 
         assert np.array_equal(raster_points, raster_points_lin)
@@ -2201,17 +2196,17 @@ def test_interp_points__synthetic(self, tag_aop: str | None, shift_aop: bool) ->
 
         points_x_in, points_y_in = raster.ij2xy(i=index_x_in, j=index_y_in, shift_area_or_point=shift_aop)
 
-        points_in = np.array((points_x_in, points_y_in)).T
-
         # Here again compare methods
-        raster_points_in = raster.interp_points(points_in, method="linear", shift_area_or_point=shift_aop)
+        raster_points_in = raster.interp_points(
+            (points_x_in, points_y_in), method="linear", shift_area_or_point=shift_aop
+        )
         raster_points_in_interpn = raster.interp_points(
-            points_in, method="linear", force_scipy_function="interpn", shift_area_or_point=shift_aop
+            (points_x_in, points_y_in), method="linear", force_scipy_function="interpn", shift_area_or_point=shift_aop
         )
 
         assert np.array_equal(raster_points_in, raster_points_in_interpn)
 
-        for i in range(len(points_in)):
+        for i in range(len(points_x_in)):
 
             xlow = int(index_x_in[i] - 0.5)
             xupp = int(index_x_in[i] + 0.5)
@@ -2228,7 +2223,8 @@ def test_interp_points__synthetic(self, tag_aop: str | None, shift_aop: bool) ->
             + [(4, i) for i in np.arange(1, 4)]
             + [(i, 4) for i in np.arange(4, 1)]
         )
-        raster_points_out = raster.interp_points(points_out)
+        points_out_xy = list(zip(*points_out))
+        raster_points_out = raster.interp_points(points_out_xy)
         assert all(~np.isfinite(raster_points_out))
 
         # To use cubic or quintic, we need a larger grid (minimum 6x6, but let's aim bigger with 50x50)
@@ -2242,14 +2238,19 @@ def test_interp_points__synthetic(self, tag_aop: str | None, shift_aop: bool) ->
         index_x_in_rand = rng.integers(low=8, high=42, size=(10,)) + rng.normal(scale=0.3)
         index_y_in_rand = rng.integers(low=8, high=42, size=(10,)) + rng.normal(scale=0.3)
         points_x_rand, points_y_rand = raster.ij2xy(i=index_x_in_rand, j=index_y_in_rand, shift_area_or_point=shift_aop)
-        points_in_rand = np.array((points_x_rand, points_y_rand)).T
 
         for method in ["nearest", "linear", "cubic", "quintic"]:
             raster_points_mapcoords = raster.interp_points(
-                points_in_rand, method=method, force_scipy_function="map_coordinates", shift_area_or_point=shift_aop
+                (points_x_rand, points_y_rand),
+                method=method,
+                force_scipy_function="map_coordinates",
+                shift_area_or_point=shift_aop,
             )
             raster_points_interpn = raster.interp_points(
-                points_in_rand, method=method, force_scipy_function="interpn", shift_area_or_point=shift_aop
+                (points_x_rand, points_y_rand),
+                method=method,
+                force_scipy_function="interpn",
+                shift_area_or_point=shift_aop,
             )
 
             # Not exactly equal in floating point precision since changes in Scipy 1.13.0,
@@ -2264,15 +2265,19 @@ def test_interp_points__synthetic(self, tag_aop: str | None, shift_aop: bool) ->
             i=index_x_edge_rand, j=index_y_edge_rand, shift_area_or_point=shift_aop
         )
 
-        points_edge_rand = np.array((points_x_rand, points_y_rand)).T
-
         # Nearest doesn't apply, just linear and above
         for method in ["cubic", "quintic"]:
             raster_points_mapcoords_edge = raster.interp_points(
-                points_edge_rand, method=method, force_scipy_function="map_coordinates", shift_area_or_point=shift_aop
+                (points_x_rand, points_y_rand),
+                method=method,
+                force_scipy_function="map_coordinates",
+                shift_area_or_point=shift_aop,
             )
             raster_points_interpn_edge = raster.interp_points(
-                points_edge_rand, method=method, force_scipy_function="interpn", shift_area_or_point=shift_aop
+                (points_x_rand, points_y_rand),
+                method=method,
+                force_scipy_function="interpn",
+                shift_area_or_point=shift_aop,
             )
 
             assert all(~np.isfinite(raster_points_mapcoords_edge))