Refactor value_at_coords into reduce_points and fix NaN behaviour (…

…#581)

doc/source/raster_class.md

@@ -302,17 +302,15 @@ prox_lt_100_from_vect = rast.proximity(vector=vect_lt_100)
 prox_lt_100_from_vect
 ```
 
-## Interpolate or extract to point
+## Interpolate or reduce to point
 
 Interpolating or extracting {class}`~geoutils.Raster` values at specific points can be done through:
-- the {func}`~geoutils.Raster.value_at_coords` function, that extracts the single closest pixel or a surrounding window for each coordinate, on which
-  can be applied reducing any function ({func}`numpy.ma.mean` by default), or
-- the {func}`~geoutils.Raster.interp_points` function, that interpolates the {class}`~geoutils.Raster`'s regular grid to each coordinate using a
-  resampling algorithm.
+- the {func}`~geoutils.Raster.reduce_points` function, that applies a reductor function ({func}`numpy.ma.mean` by default) to a surrounding window for each coordinate, or
+- the {func}`~geoutils.Raster.interp_points` function, that interpolates the {class}`~geoutils.Raster`'s regular grid to each coordinate using a resampling algorithm.
 
 ```{code-cell} ipython3
 # Extract median value in a 3 x 3 pixel window
-rast_reproj.value_at_coords(x=0.5, y=0.5, window=3, reducer_function=np.ma.median)
+rast_reproj.reduce_points((0.5, 0.5), window=3, reducer_function=np.ma.median)
 ```
 
 ```{code-cell} ipython3
@@ -321,7 +319,7 @@ rast_reproj.interp_points((0.5, 0.5), method="quintic")
 ```
 
 ```{note}
-Both {func}`~geoutils.Raster.value_at_coords` and {func}`~geoutils.Raster.interp_points` can be passed a single coordinate as {class}`floats<float>`, or a
+Both {func}`~geoutils.Raster.reduce_points` and {func}`~geoutils.Raster.interp_points` can be passed a single coordinate as {class}`floats<float>`, or a
 {class}`list` of coordinates.
 ```
 

examples/analysis/point_extraction/reduction.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.value_at_coords(x=x_coords, y=y_coords)
+vals = rast.reduce_points((x_coords, y_coords))
 
 # %%
 # Replace by Vector function once done
@@ -40,7 +40,7 @@
 # By default, :func:`~geoutils.Raster.value_at_coords` extracts the closest pixel value. But it can also be passed a window size and reductor function to
 # extract an average value or other statistic based on neighbouring pixels.
 
-vals_reduced = rast.value_at_coords(x=x_coords, y=y_coords, window=5, reducer_function=np.nanmedian)
+vals_reduced = rast.reduce_points((x_coords, y_coords), window=5, reducer_function=np.nanmedian)
 
 np.nanmean(vals - vals_reduced)
 
@@ -50,8 +50,8 @@
 
 # Replace by Vector fonction once done
 coords = rast.coords(grid=True)
-x_closest = rast.copy(new_array=coords[0]).value_at_coords(x=x_coords, y=y_coords).squeeze()
-y_closest = rast.copy(new_array=coords[1]).value_at_coords(x=x_coords, y=y_coords).squeeze()
+x_closest = rast.copy(new_array=coords[0]).reduce_points((x_coords, y_coords)).squeeze()
+y_closest = rast.copy(new_array=coords[1]).reduce_points((x_coords, y_coords)).squeeze()
 from shapely.geometry import box
 
 geometry = [

geoutils/raster/raster.py

@@ -49,6 +49,7 @@
     _get_bounds_projected,
     _get_footprint_projected,
     _get_utm_ups_crs,
+    reproject_from_latlon,
 )
 from geoutils.raster.array import get_mask_from_array
 from geoutils.raster.sampling import subsample_array
@@ -3301,31 +3302,32 @@ def plot(
         else:
             return None
 
-    def value_at_coords(
+    def reduce_points(
         self,
-        x: Number | ArrayLike,
-        y: Number | ArrayLike,
-        latlon: bool = False,
+        points: tuple[ArrayLike, ArrayLike],
+        reducer_function: Callable[[NDArrayNum], float] = np.ma.mean,
+        window: int | None = None,
+        input_latlon: bool = False,
         band: int | None = None,
         masked: bool = False,
-        window: int | None = None,
-        reducer_function: Callable[[NDArrayNum], float] = np.ma.mean,
         return_window: bool = False,
         boundless: bool = True,
     ) -> Any:
         """
-        Extract raster values at the nearest pixels from the specified coordinates,
-        or reduced (e.g., mean of pixels) from a window around the specified coordinates.
+        Reduce raster values around point coordinates.
+
 
         By default, samples pixel value of each band. Can be passed a band index to sample from.
 
-        :param x: X (or longitude) coordinate(s).
-        :param y: Y (or latitude) coordinate(s).
-        :param latlon: Whether coordinates are provided as longitude-latitude.
+        Uses Rasterio's windowed reading to keep memory usage low (for a raster not loaded).
+
+        :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 reducer_function: Reducer function to apply to the values in window (defaults to np.mean).
+        :param window: Window size to read around coordinates. Must be odd.
+        :param input_latlon: Whether the input is in latlon, unregarding of Raster CRS
         :param band: Band number to extract from (from 1 to self.count).
         :param masked: Whether to return a masked array, or classic array.
-        :param window: Window size to read around coordinates. Must be odd.
-        :param reducer_function: Reducer function to apply to the values in window (defaults to np.mean).
         :param return_window: Whether to return the windows (in addition to the reduced value).
         :param boundless: Whether to allow windows that extend beyond the extent.
 
@@ -3344,6 +3346,9 @@ def value_at_coords(
             (c = provided coordinate, v= value of surrounding coordinate)
 
         """
+
+        x, y = points
+
         # Check for array-like inputs
         if (
             not isinstance(x, (float, np.floating, int, np.integer))
@@ -3391,10 +3396,8 @@ def format_value(value: Any) -> Any:
             list_windows = []
 
         # Convert to latlon if asked
-        if latlon:
-            from geoutils import projtools
-
-            x, y = projtools.reproject_from_latlon((y, x), self.crs)  # type: ignore
+        if input_latlon:
+            x, y = reproject_from_latlon((y, x), self.crs)  # type: ignore
 
         # Convert coordinates to pixel space
         rows, cols = rio.transform.rowcol(self.transform, x, y, op=floor)
@@ -3433,7 +3436,7 @@ def format_value(value: Any) -> Any:
                 else:
                     data = self.data[slice(None) if band is None else band - 1, row : row + height, col : col + width]
                 if not masked:
-                    data = data.filled()
+                    data = data.astype(np.float32).filled(np.nan)
                 value = format_value(data)
                 win: NDArrayNum | dict[int, NDArrayNum] = data
 
@@ -3685,7 +3688,7 @@ def interp_points(
         :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).
-        :param input_latlon: Whether the input is in latlon, unregarding of Raster CRS
+        :param input_latlon: Whether the input is in latlon, unregarding of Raster CRS.
         :param shift_area_or_point: Whether to shift with pixel interpretation, which shifts to center of pixel
             coordinates if self.area_or_point is "Point" and maintains corner pixel coordinate if it is "Area" or None.
             Defaults to True. Can be configured with the global setting geoutils.config["shift_area_or_point"].

tests/test_raster/test_raster.py

@@ -2295,12 +2295,12 @@ def test_interp_points__real(self) -> None:
         val_latlon = r.interp_points((lat, lon), method="linear", input_latlon=True)[0]
         assert val == pytest.approx(val_latlon, abs=0.0001)
 
-    def test_value_at_coords(self) -> None:
+    def test_reduce_points(self) -> None:
         """
-        Test that value at coords works as intended
+        Test reduce points.
         """
 
-        # -- Tests 1: check based on indexed values --
+        # -- Tests 1: Check based on indexed values --
 
         # Open raster
         r = gu.Raster(self.landsat_b4_crop_path)
@@ -2321,21 +2321,21 @@ def test_value_at_coords(self) -> None:
         assert y_out == ytest0
 
         # Check that the value at this coordinate is the same as when indexing
-        z_val = r.value_at_coords(xtest0, ytest0)
+        z_val = r.reduce_points((xtest0, ytest0))
         z = r.data.data[itest0, jtest0]
         assert z == z_val
 
         # Check that the value is the same the other 4 corners of the pixel
-        assert z == r.value_at_coords(xtest0 + 0.49 * r.res[0], ytest0 - 0.49 * r.res[1])
-        assert z == r.value_at_coords(xtest0 - 0.49 * r.res[0], ytest0 + 0.49 * r.res[1])
-        assert z == r.value_at_coords(xtest0 - 0.49 * r.res[0], ytest0 - 0.49 * r.res[1])
-        assert z == r.value_at_coords(xtest0 + 0.49 * r.res[0], ytest0 + 0.49 * r.res[1])
+        assert z == r.reduce_points((xtest0 + 0.49 * r.res[0], ytest0 - 0.49 * r.res[1]))
+        assert z == r.reduce_points((xtest0 - 0.49 * r.res[0], ytest0 + 0.49 * r.res[1]))
+        assert z == r.reduce_points((xtest0 - 0.49 * r.res[0], ytest0 - 0.49 * r.res[1]))
+        assert z == r.reduce_points((xtest0 + 0.49 * r.res[0], ytest0 + 0.49 * r.res[1]))
 
         # -- Tests 2: check arguments work as intended --
 
         # 1/ Lat-lon argument check by getting the coordinates of our last test point
         lat, lon = reproject_to_latlon(points=[xtest0, ytest0], in_crs=r.crs)
-        z_val_2 = r.value_at_coords(lon, lat, latlon=True)
+        z_val_2 = r.reduce_points((lon, lat), input_latlon=True)
         assert z_val == z_val_2
 
         # 2/ Band argument
@@ -2345,22 +2345,22 @@ def test_value_at_coords(self) -> None:
         itest = int(itest[0])
         jtest = int(jtest[0])
         # Extract the values
-        z_band1 = r_multi.value_at_coords(xtest0, ytest0, band=1)
-        z_band2 = r_multi.value_at_coords(xtest0, ytest0, band=2)
-        z_band3 = r_multi.value_at_coords(xtest0, ytest0, band=3)
+        z_band1 = r_multi.reduce_points((xtest0, ytest0), band=1)
+        z_band2 = r_multi.reduce_points((xtest0, ytest0), band=2)
+        z_band3 = r_multi.reduce_points((xtest0, ytest0), band=3)
         # Compare to the Raster array slice
         assert list(r_multi.data[:, itest, jtest]) == [z_band1, z_band2, z_band3]
 
         # 3/ Masked argument
         r_multi.data[:, itest, jtest] = np.ma.masked
-        z_not_ma = r_multi.value_at_coords(xtest0, ytest0, band=1)
+        z_not_ma = r_multi.reduce_points((xtest0, ytest0), band=1)
         assert not np.ma.is_masked(z_not_ma)
-        z_ma = r_multi.value_at_coords(xtest0, ytest0, band=1, masked=True)
+        z_ma = r_multi.reduce_points((xtest0, ytest0), band=1, masked=True)
         assert np.ma.is_masked(z_ma)
 
         # 4/ Window argument
-        val_window, z_window = r_multi.value_at_coords(
-            xtest0, ytest0, band=1, window=3, masked=True, return_window=True
+        val_window, z_window = r_multi.reduce_points(
+            (xtest0, ytest0), band=1, window=3, masked=True, return_window=True
         )
         assert (
             val_window
@@ -2370,37 +2370,37 @@ def test_value_at_coords(self) -> None:
         assert np.array_equal(z_window, r_multi.data[0, itest - 1 : itest + 2, jtest - 1 : jtest + 2])
 
         # 5/ Reducer function argument
-        val_window2 = r_multi.value_at_coords(
-            xtest0, ytest0, band=1, window=3, masked=True, reducer_function=np.ma.median
+        val_window2 = r_multi.reduce_points(
+            (xtest0, ytest0), band=1, window=3, masked=True, reducer_function=np.ma.median
         )
         assert val_window2 == np.ma.median(r_multi.data[0, itest - 1 : itest + 2, jtest - 1 : jtest + 2])
 
         # -- Tests 3: check that errors are raised when supposed for non-boolean arguments --
 
         # Verify that passing a window that is not a whole number fails
         with pytest.raises(ValueError, match=re.escape("Window must be a whole number.")):
-            r.value_at_coords(xtest0, ytest0, window=3.5)  # type: ignore
+            r.reduce_points((xtest0, ytest0), window=3.5)  # type: ignore
         # Same for an odd number
         with pytest.raises(ValueError, match=re.escape("Window must be an odd number.")):
-            r.value_at_coords(xtest0, ytest0, window=4)
+            r.reduce_points((xtest0, ytest0), window=4)
         # But a window that is a whole number as a float works
-        r.value_at_coords(xtest0, ytest0, window=3.0)  # type: ignore
+        r.reduce_points((xtest0, ytest0), window=3.0)  # type: ignore
 
         # -- Tests 4: check that passing an array-like object works
 
         # For simple coordinates
         x_coords = [xtest0, xtest0 + 100]
         y_coords = [ytest0, ytest0 - 100]
-        vals = r_multi.value_at_coords(x=x_coords, y=y_coords)
-        val0, win0 = r_multi.value_at_coords(x=x_coords[0], y=y_coords[0], return_window=True)
-        val1, win1 = r_multi.value_at_coords(x=x_coords[1], y=y_coords[1], return_window=True)
+        vals = r_multi.reduce_points((x_coords, y_coords))
+        val0, win0 = r_multi.reduce_points((x_coords[0], y_coords[0]), return_window=True)
+        val1, win1 = r_multi.reduce_points((x_coords[1], y_coords[1]), return_window=True)
 
         assert len(vals) == len(x_coords)
-        assert vals[0] == val0
-        assert vals[1] == val1
+        assert np.array_equal(vals[0], val0, equal_nan=True)
+        assert np.array_equal(vals[1], val1, equal_nan=True)
 
         # With a return window argument
-        vals, windows = r_multi.value_at_coords(x=x_coords, y=y_coords, return_window=True)
+        vals, windows = r_multi.reduce_points((x_coords, y_coords), return_window=True)
         assert len(windows) == len(x_coords)
         assert np.array_equal(windows[0], win0, equal_nan=True)
         assert np.array_equal(windows[1], win1, equal_nan=True)
@@ -2410,17 +2410,17 @@ def test_value_at_coords(self) -> None:
         # Lower right pixel
         x, y = [r.bounds.right - r.res[0] / 2, r.bounds.bottom + r.res[1] / 2]
         lat, lon = pt.reproject_to_latlon([x, y], r.crs)
-        assert r.value_at_coords(x, y) == r.value_at_coords(lon, lat, latlon=True) == r.data[-1, -1]
+        assert r.reduce_points((x, y)) == r.reduce_points((lon, lat), input_latlon=True) == r.data[-1, -1]
 
         # One pixel above
         x, y = [r.bounds.right - r.res[0] / 2, r.bounds.bottom + 3 * r.res[1] / 2]
         lat, lon = pt.reproject_to_latlon([x, y], r.crs)
-        assert r.value_at_coords(x, y) == r.value_at_coords(lon, lat, latlon=True) == r.data[-2, -1]
+        assert r.reduce_points((x, y)) == r.reduce_points((lon, lat), input_latlon=True) == r.data[-2, -1]
 
         # One pixel left
         x, y = [r.bounds.right - 3 * r.res[0] / 2, r.bounds.bottom + r.res[1] / 2]
         lat, lon = pt.reproject_to_latlon([x, y], r.crs)
-        assert r.value_at_coords(x, y) == r.value_at_coords(lon, lat, latlon=True) == r.data[-1, -2]
+        assert r.reduce_points((x, y)) == r.reduce_points((lon, lat), input_latlon=True) == r.data[-1, -2]
 
     @pytest.mark.parametrize("example", [landsat_b4_path, aster_dem_path])  # type: ignore
     def test_set_nodata(self, example: str) -> None: