Skip to content

Commit

Permalink
Add bounds to georeferencing, accept None input for interpolate with …
Browse files Browse the repository at this point in the history
…return_interpolator
  • Loading branch information
rhugonnet committed Aug 23, 2024
1 parent d36b713 commit de39c35
Show file tree
Hide file tree
Showing 3 changed files with 19 additions and 11 deletions.
6 changes: 6 additions & 0 deletions geoutils/raster/georeferencing.py
Original file line number Diff line number Diff line change
Expand Up @@ -164,3 +164,9 @@ def _res(transform: rio.transform.Affine) -> tuple[float, float]:
"""See description of Raster.res"""

return transform[0], abs(transform[4])


def _bounds(transform: rio.transform.Affine, shape: tuple[int, int]):
"""See description of Raster.bounds."""

return rio.coords.BoundingBox(*rio.transform.array_bounds(height=shape[0], width=shape[1], transform=transform))
20 changes: 11 additions & 9 deletions geoutils/raster/interpolate.py
Original file line number Diff line number Diff line change
Expand Up @@ -237,7 +237,7 @@ def _interp_points(
array: NDArrayNum,
transform: rio.transform.Affine,
area_or_point: Literal["Area", "Point"] | None,
points: tuple[Number, Number] | tuple[NDArrayNum, NDArrayNum],
points: tuple[Number, Number] | tuple[NDArrayNum, NDArrayNum] | None,
method: Literal["nearest", "linear", "cubic", "quintic", "slinear", "pchip", "splinef2d"] = "linear",
dist_nodata_spread: Literal["half_order_up", "half_order_down"] | int = "half_order_up",
shift_area_or_point: bool | None = None,
Expand All @@ -250,16 +250,18 @@ def _interp_points(
# If array is not a floating dtype (to support NaNs), convert dtype
if not np.issubdtype(array.dtype, np.floating):
array = array.astype(np.float32)
shape: tuple[int, int] = array.shape[0:2] # type: ignore

# Get coordinates
x, y = points

i, j = _xy2ij(x, y, transform=transform, area_or_point=area_or_point, shift_area_or_point=shift_area_or_point)
# TODO: Add check about None for "points" depending on "return_interpolator"
if not return_interpolator:
x, y = points
i, j = _xy2ij(x, y, transform=transform, area_or_point=area_or_point,
shift_area_or_point=shift_area_or_point)

shape: tuple[int, int] = array.shape[0:2] # type: ignore
ind_invalid = np.vectorize(
lambda k1, k2: _outside_image(k1, k2, transform=transform, area_or_point=area_or_point, shape=shape, index=True)
)(j, i)
ind_invalid = np.vectorize(
lambda k1, k2: _outside_image(k1, k2, transform=transform, area_or_point=area_or_point, shape=shape,
index=True)
)(j, i)

# If the raster is on an equal grid, use scipy.ndimage.map_coordinates
force_map_coords = force_scipy_function is not None and force_scipy_function == "map_coordinates"
Expand Down
4 changes: 2 additions & 2 deletions geoutils/raster/raster.py
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@
reproject_from_latlon,
)
from geoutils.raster.array import get_mask_from_array
from geoutils.raster.georeferencing import _coords, _ij2xy, _outside_image, _res, _xy2ij
from geoutils.raster.georeferencing import _coords, _ij2xy, _outside_image, _res, _xy2ij, _bounds
from geoutils.raster.interpolate import _interp_points
from geoutils.raster.sampling import subsample_array
from geoutils.vector import Vector
Expand Down Expand Up @@ -793,7 +793,7 @@ def res(self) -> tuple[float | int, float | int]:
@property
def bounds(self) -> rio.coords.BoundingBox:
"""Bounding coordinates of the raster."""
return rio.coords.BoundingBox(*rio.transform.array_bounds(self.height, self.width, self.transform))
return _bounds(transform=self.transform, shape=self.shape)

@property
def footprint(self) -> Vector:
Expand Down

0 comments on commit de39c35

Please sign in to comment.