From dd7e6f71ed9d11635b175824aafb980153d53491 Mon Sep 17 00:00:00 2001 From: Romain Hugonnet Date: Thu, 30 May 2024 15:56:31 -0800 Subject: [PATCH] Small fixes and optimizations for `to_pointcloud` and `interp_points` (#549) --- doc/source/raster_class.md | 2 +- .../point_extraction/interpolation.py | 4 +- geoutils/raster/raster.py | 17 ++++-- tests/test_raster/test_raster.py | 59 ++++++++++--------- 4 files changed, 46 insertions(+), 36 deletions(-) diff --git a/doc/source/raster_class.md b/doc/source/raster_class.md index f31c22a4..379c33cc 100644 --- a/doc/source/raster_class.md +++ b/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} diff --git a/examples/analysis/point_extraction/interpolation.py b/examples/analysis/point_extraction/interpolation.py index bef62805..0a89cec9 100644 --- a/examples/analysis/point_extraction/interpolation.py +++ b/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) # %% diff --git a/geoutils/raster/raster.py b/geoutils/raster/raster.py index bc05d62b..c127c75c 100644 --- a/geoutils/raster/raster.py +++ b/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 diff --git a/tests/test_raster/test_raster.py b/tests/test_raster/test_raster.py index fa32c440..c4bd1f17 100644 --- a/tests/test_raster/test_raster.py +++ b/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,7 +2131,7 @@ 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) @@ -2142,13 +2139,13 @@ def test_xy2ij(self) -> None: 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))