Skip to content

Commit

Permalink
Fix interpn option of interp_points (#554)
Browse files Browse the repository at this point in the history
  • Loading branch information
rhugonnet authored Jun 1, 2024
1 parent dd7e6f7 commit 86bc449
Show file tree
Hide file tree
Showing 2 changed files with 28 additions and 15 deletions.
7 changes: 5 additions & 2 deletions geoutils/raster/raster.py
Original file line number Diff line number Diff line change
Expand Up @@ -3729,7 +3729,7 @@ def interp_points(
# Otherwise, use scipy.interpolate.interpn
else:
# Get lower-left corner coordinates
xycoords = self.coords(grid=False, force_offset="ll")
xycoords = self.coords(grid=False, shift_area_or_point=shift_area_or_point)

# Let interpolation outside the bounds not raise any error by default
if "bounds_error" not in kwargs.keys():
Expand All @@ -3738,7 +3738,10 @@ def interp_points(
if "fill_value" not in kwargs.keys():
kwargs.update({"fill_value": np.nan})

rpoints = interpn(xycoords, self.get_nanarray(), np.array([i, j]).T, method=method, **kwargs)
# Using direct coordinates, Y is the first axis, and we need to flip it
rpoints = interpn(
(np.flip(xycoords[1], axis=0), xycoords[0]), self.get_nanarray(), (y, x), method=method, **kwargs
)

rpoints = np.array(rpoints, dtype=np.float32)
rpoints[np.array(ind_invalid)] = np.nan
Expand Down
36 changes: 23 additions & 13 deletions tests/test_raster/test_raster.py
Original file line number Diff line number Diff line change
Expand Up @@ -2135,19 +2135,6 @@ def test_xy2ij(self) -> None:

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_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]
assert val == pytest.approx(val_latlon, abs=0.0001)

@pytest.mark.parametrize("tag_aop", [None, "Area", "Point"]) # type: ignore
@pytest.mark.parametrize("shift_aop", [True, False]) # type: ignore
def test_interp_points__synthetic(self, tag_aop: str | None, shift_aop: bool) -> None:
Expand Down Expand Up @@ -2283,6 +2270,29 @@ def test_interp_points__synthetic(self, tag_aop: str | None, shift_aop: bool) ->
assert all(~np.isfinite(raster_points_mapcoords_edge))
assert all(~np.isfinite(raster_points_interpn_edge))

def test_interp_points__real(self) -> None:
"""Test interp_points for real data."""

r = gu.Raster(self.landsat_b4_path)
r.set_area_or_point("Area", shift_area_or_point=False)

# 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_img = r.data[int(i[0]), int(j[0])]
assert val_img == val

# Check the result is exactly the same for both methods
val2 = r.interp_points((x, y), method="linear", force_scipy_function="interpn")[0]
assert val2 == pytest.approx(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]
assert val == pytest.approx(val_latlon, abs=0.0001)

def test_value_at_coords(self) -> None:
"""
Test that value at coords works as intended
Expand Down

0 comments on commit 86bc449

Please sign in to comment.