diff --git a/geoutils/raster/raster.py b/geoutils/raster/raster.py
index c127c75c..b3661ea4 100644
--- a/geoutils/raster/raster.py
+++ b/geoutils/raster/raster.py
@@ -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():
@@ -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
diff --git a/tests/test_raster/test_raster.py b/tests/test_raster/test_raster.py
index c4bd1f17..58b54f18 100644
--- a/tests/test_raster/test_raster.py
+++ b/tests/test_raster/test_raster.py
@@ -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:
@@ -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