diff --git a/doc/source/raster_class.md b/doc/source/raster_class.md index 379c33cc..ae8b6475 100644 --- a/doc/source/raster_class.md +++ b/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`, or a +Both {func}`~geoutils.Raster.reduce_points` and {func}`~geoutils.Raster.interp_points` can be passed a single coordinate as {class}`floats`, or a {class}`list` of coordinates. ``` diff --git a/examples/analysis/point_extraction/reduction.py b/examples/analysis/point_extraction/reduction.py index b76e2c64..d9f68d3e 100644 --- a/examples/analysis/point_extraction/reduction.py +++ b/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 = [ diff --git a/geoutils/raster/raster.py b/geoutils/raster/raster.py index 43f1fe9d..c817c6dd 100644 --- a/geoutils/raster/raster.py +++ b/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"]. diff --git a/tests/test_raster/test_raster.py b/tests/test_raster/test_raster.py index c56c27d7..da6334b2 100644 --- a/tests/test_raster/test_raster.py +++ b/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,8 +2370,8 @@ 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]) @@ -2379,28 +2379,28 @@ def test_value_at_coords(self) -> None: # 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: