diff --git a/README.md b/README.md
index f5ca0143..79292f6f 100644
--- a/README.md
+++ b/README.md
@@ -77,7 +77,7 @@ Please see our [contribution page](CONTRIBUTING.md) for more detailed instructio
 See the documentation at https://xdem.readthedocs.io
 
 ### Testing - again please read!
-These tools are only valuable if we can rely on them to perform exactly as we expect. So, we need testing. Please create tests for every function that you make, as much as you are able. Guidance/examples here for the moment: https://github.com/GeoUtils/georaster/blob/master/test/test_georaster.py
+These tools are only valuable if we can rely on them to perform exactly as we expect. So, we need testing. Please create tests for every function that you make, as much as you are able. Guidance/examples here for the moment: https://github.com/GeoUtils/raster/blob/master/test/test_raster.py
 https://github.com/corteva/rioxarray/blob/master/test/integration/test_integration__io.py
 
 
diff --git a/dev-environment.yml b/dev-environment.yml
index b54fa046..01ec65eb 100644
--- a/dev-environment.yml
+++ b/dev-environment.yml
@@ -37,5 +37,5 @@ dependencies:
   - richdem
 
   - pip:
-    # - git+https://github.com/GlacioHack/GeoUtils.git
-    - geoutils==0.0.10
+     - git+https://github.com/GlacioHack/GeoUtils.git
+#    - geoutils==0.0.10
diff --git a/docs/source/code/spatialstats.py b/docs/source/code/spatialstats.py
index c012a9ac..f1774cb0 100644
--- a/docs/source/code/spatialstats.py
+++ b/docs/source/code/spatialstats.py
@@ -5,20 +5,20 @@
 import xdem
 
 # Load data
-dh = gu.georaster.Raster(xdem.examples.get_path("longyearbyen_ddem"))
+dh = gu.Raster(xdem.examples.get_path("longyearbyen_ddem"))
 ref_dem = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"))
-glacier_mask = gu.geovector.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines"))
+glacier_mask = gu.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines"))
 mask = glacier_mask.create_mask(dh)
 
 slope = xdem.terrain.get_terrain_attribute(ref_dem, attribute=["slope"])
 
 # Keep only stable terrain data
 dh.set_mask(mask)
-dh_arr = gu.spatial_tools.get_array_and_mask(dh)[0]
-slope_arr = gu.spatial_tools.get_array_and_mask(slope)[0]
+dh_arr = gu.raster.get_array_and_mask(dh)[0]
+slope_arr = gu.raster.get_array_and_mask(slope)[0]
 
 # Subsample to run the snipped code faster
-indices = gu.spatial_tools.subsample_raster(dh_arr, subsample=10000, return_indices=True, random_state=42)
+indices = gu.raster.subsample_array(dh_arr, subsample=10000, return_indices=True, random_state=42)
 dh_arr = dh_arr[indices]
 slope_arr = slope_arr[indices]
 
diff --git a/docs/source/code/spatialstats_heterosc_slope.py b/docs/source/code/spatialstats_heterosc_slope.py
index a2cd98ca..d74a211f 100644
--- a/docs/source/code/spatialstats_heterosc_slope.py
+++ b/docs/source/code/spatialstats_heterosc_slope.py
@@ -5,9 +5,9 @@
 import xdem
 
 # Load data
-dh = gu.georaster.Raster(xdem.examples.get_path("longyearbyen_ddem"))
+dh = gu.Raster(xdem.examples.get_path("longyearbyen_ddem"))
 ref_dem = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"))
-glacier_mask = gu.geovector.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines"))
+glacier_mask = gu.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines"))
 mask = glacier_mask.create_mask(dh)
 
 # Get slope for non-stationarity
diff --git a/environment.yml b/environment.yml
index 72e3f30f..e71c9620 100644
--- a/environment.yml
+++ b/environment.yml
@@ -23,5 +23,5 @@ dependencies:
   - pip
 
   - pip:
-    # - git+https://github.com/GlacioHack/GeoUtils.git
-    - geoutils==0.0.10
+     - git+https://github.com/GlacioHack/GeoUtils.git
+#    - geoutils==0.0.10
diff --git a/examples/advanced/plot_heterosc_estimation_modelling.py b/examples/advanced/plot_heterosc_estimation_modelling.py
index 0f72cc16..f94a825d 100644
--- a/examples/advanced/plot_heterosc_estimation_modelling.py
+++ b/examples/advanced/plot_heterosc_estimation_modelling.py
@@ -44,11 +44,11 @@
 
 # %%
 # We convert to arrays and keep only stable terrain for the analysis of variability
-dh_arr = gu.spatial_tools.get_array_and_mask(dh)[0][~mask_glacier]
-slope_arr = gu.spatial_tools.get_array_and_mask(slope)[0][~mask_glacier]
-aspect_arr = gu.spatial_tools.get_array_and_mask(aspect)[0][~mask_glacier]
-planc_arr = gu.spatial_tools.get_array_and_mask(planc)[0][~mask_glacier]
-profc_arr = gu.spatial_tools.get_array_and_mask(profc)[0][~mask_glacier]
+dh_arr = dh.get_nanarray()[~mask_glacier]
+slope_arr = slope.get_nanarray()[~mask_glacier]
+aspect_arr = aspect.get_nanarray()[~mask_glacier]
+planc_arr = planc.get_nanarray()[~mask_glacier]
+profc_arr = profc.get_nanarray()[~mask_glacier]
 
 # %%
 # We use :func:`xdem.spatialstats.nd_binning` to perform N-dimensional binning on all those terrain variables, with uniform
diff --git a/examples/basic/plot_dem_subtraction.py b/examples/basic/plot_dem_subtraction.py
index 8b688100..289c547f 100644
--- a/examples/basic/plot_dem_subtraction.py
+++ b/examples/basic/plot_dem_subtraction.py
@@ -38,7 +38,7 @@
 # To hide this prompt, add ``.reproject(..., silent=True)``.
 # By default, :func:`xdem.DEM.reproject` uses "bilinear" resampling (assuming resampling is needed).
 # Other options are "nearest" (fast but inaccurate), "cubic_spline", "lanczos" and others.
-# See `geoutils.Raster.reproject() <https://geoutils.readthedocs.io/en/latest/api.html#geoutils.georaster.Raster.reproject>`_ and `rasterio.enums.Resampling <https://rasterio.readthedocs.io/en/latest/api/rasterio.enums.html#rasterio.enums.Resampling>`_ for more information about reprojection.
+# See `geoutils.Raster.reproject() <https://geoutils.readthedocs.io/en/latest/api.html#geoutils.raster.Raster.reproject>`_ and `rasterio.enums.Resampling <https://rasterio.readthedocs.io/en/latest/api/rasterio.enums.html#rasterio.enums.Resampling>`_ for more information about reprojection.
 #
 # Now, we are ready to generate the dDEM:
 
diff --git a/tests/test_coreg.py b/tests/test_coreg.py
index d09556af..79e25133 100644
--- a/tests/test_coreg.py
+++ b/tests/test_coreg.py
@@ -14,7 +14,7 @@
 import pytest
 import rasterio as rio
 from geoutils import Raster, Vector
-from geoutils.georaster.raster import RasterType
+from geoutils.raster import RasterType
 
 with warnings.catch_warnings():
     warnings.simplefilter("ignore")
@@ -264,12 +264,12 @@ def test_deramping(self) -> None:
         deramp.fit(**self.fit_params)
 
         # Apply the deramping to a DEM
-        deramped_dem, _ = deramp.apply(self.tba.data, self.ref.transform, self.ref.crs)
+        deramped_dem = deramp.apply(self.tba)
 
         # Get the periglacial offset after deramping
-        periglacial_offset = (self.ref.data.squeeze() - deramped_dem)[self.inlier_mask.squeeze()]
+        periglacial_offset = (self.ref - deramped_dem)[self.inlier_mask]
         # Get the periglacial offset before deramping
-        pre_offset = ((self.ref.data - self.tba.data)[self.inlier_mask]).squeeze()
+        pre_offset = (self.ref - self.tba)[self.inlier_mask]
 
         # Check that the error improved
         assert np.abs(np.mean(periglacial_offset)) < np.abs(np.mean(pre_offset))
@@ -434,7 +434,7 @@ def test_z_scale_corr(self) -> None:
         assert np.abs(np.nanmedian(diff)) < 0.01
 
         # Create a spatially correlated error field to mess with the algorithm a bit.
-        corr_size = int(self.ref.data.shape[2] / 100)
+        corr_size = int(self.ref.data.shape[1] / 100)
         error_field = cv2.resize(
             cv2.GaussianBlur(
                 np.repeat(
@@ -442,7 +442,7 @@ def test_z_scale_corr(self) -> None:
                         np.random.randint(
                             0,
                             255,
-                            (self.ref.data.shape[1] // corr_size, self.ref.data.shape[2] // corr_size),
+                            (self.ref.data.shape[0] // corr_size, self.ref.data.shape[1] // corr_size),
                             dtype="uint8",
                         ),
                         corr_size,
@@ -455,7 +455,7 @@ def test_z_scale_corr(self) -> None:
                 sigmaX=corr_size,
             )
             / 255,
-            dsize=(self.ref.data.shape[2], self.ref.data.shape[1]),
+            dsize=(self.ref.data.shape[1], self.ref.data.shape[0]),
         )
 
         # Create 50000 random nans
@@ -517,13 +517,13 @@ def test_blockwise_coreg(self, pipeline: coreg.Coreg, subdivision: int) -> None:
         chunk_numbers = [m["i"] for m in blockwise._meta["coreg_meta"]]
         assert np.unique(chunk_numbers).shape[0] == len(chunk_numbers)
 
-        transformed_dem, _ = blockwise.apply(self.tba.data, self.tba.transform, self.tba.crs)
+        transformed_dem = blockwise.apply(self.tba)
 
-        ddem_pre = (self.ref.data - self.tba.data)[~self.inlier_mask].squeeze().filled(np.nan)
-        ddem_post = (self.ref.data.squeeze() - transformed_dem)[~self.inlier_mask.squeeze()].filled(np.nan)
+        ddem_pre = (self.ref - self.tba)[~self.inlier_mask]
+        ddem_post = (self.ref - transformed_dem)[~self.inlier_mask]
 
         # Check that the periglacial difference is lower after coregistration.
-        assert abs(np.nanmedian(ddem_post)) < abs(np.nanmedian(ddem_pre))
+        assert abs(np.ma.median(ddem_post)) < abs(np.ma.median(ddem_pre))
 
         stats = blockwise.stats()
 
@@ -554,7 +554,7 @@ def test_blockwise_coreg_large_gaps(self) -> None:
         # Copy the TBA DEM and set a square portion to nodata
         tba = self.tba.copy()
         mask = np.zeros(np.shape(tba.data), dtype=bool)
-        mask[0, 450:500, 450:500] = True
+        mask[450:500, 450:500] = True
         tba.set_mask(mask=mask)
 
         blockwise = xdem.coreg.BlockwiseCoreg(xdem.coreg.NuthKaab(), 8, warn_failures=False)
@@ -583,7 +583,7 @@ def test_coreg_raster_and_ndarray_args(self) -> None:
             nodata=-9999,
         )
         # Assign a funny value to one particular pixel. This is to validate that reprojection works perfectly.
-        dem1.data[0, 1, 1] = 100
+        dem1.data[1, 1] = 100
 
         # Translate the DEM 1 "meter" right and add a bias
         dem2 = dem1.reproject(dst_bounds=rio.coords.BoundingBox(1, 0, 6, 5), silent=True)
@@ -822,7 +822,7 @@ def test_coreg_oneliner(self) -> None:
 def test_apply_matrix() -> None:
     warnings.simplefilter("error")
     ref, tba, outlines = load_examples()  # Load example reference, to-be-aligned and mask.
-    ref_arr = gu.spatial_tools.get_array_and_mask(ref)[0]
+    ref_arr = gu.raster.get_array_and_mask(ref)[0]
 
     # Test only bias (it should just apply the bias and not make anything else)
     bias = 5
@@ -1043,7 +1043,7 @@ def test_create_inlier_mask() -> None:
 
     # - Assert that without filtering create_inlier_mask behaves as if calling Vector.create_mask - #
     # Masking inside - using Vector
-    inlier_mask_comp = ~outlines.create_mask(ref)
+    inlier_mask_comp = ~outlines.create_mask(ref, as_array=True)
     inlier_mask = xdem.coreg.create_inlier_mask(
         tba,
         ref,
@@ -1253,7 +1253,7 @@ def test_dem_coregistration() -> None:
         ],
         resample=True,
     )
-    gl_mask = outlines.create_mask(dem_coreg)
+    gl_mask = outlines.create_mask(dem_coreg, as_array=True)
     assert np.all(~inlier_mask[gl_mask])
 
     # Testing with plot
diff --git a/tests/test_ddem.py b/tests/test_ddem.py
index d9f80906..1bddc812 100644
--- a/tests/test_ddem.py
+++ b/tests/test_ddem.py
@@ -29,7 +29,7 @@ def test_copy(self) -> None:
 
         ddem2.data += 1
 
-        assert self.ddem != ddem2
+        assert not self.ddem.raster_equal(ddem2)
 
     def test_filled_data(self) -> None:
         """Test that the filled_data property points to the right data."""
diff --git a/tests/test_dem.py b/tests/test_dem.py
index 3a3890db..aaf5b301 100644
--- a/tests/test_dem.py
+++ b/tests/test_dem.py
@@ -2,13 +2,13 @@
 import os
 import warnings
 
-import geoutils.georaster as gr
-import geoutils.satimg as si
+import geoutils.raster as gr
+import geoutils.raster.satimg as si
 import numpy as np
 import pyproj
 import pytest
 import rasterio as rio
-from geoutils.georaster.raster import _default_rio_attrs
+from geoutils.raster.raster import _default_rio_attrs
 
 import xdem
 from xdem.dem import DEM
@@ -94,7 +94,7 @@ def test_copy(self) -> None:
         assert isinstance(r2, xdem.dem.DEM)
 
         # Check all immutable attributes are equal
-        # georaster_attrs = ['bounds', 'count', 'crs', 'dtypes', 'height', 'indexes', 'nodata',
+        # raster_attrs = ['bounds', 'count', 'crs', 'dtypes', 'height', 'indexes', 'nodata',
         #                    'res', 'shape', 'transform', 'width']
         # satimg_attrs = ['satellite', 'sensor', 'product', 'version', 'tile_name', 'datetime']
         # dem_attrs = ['vref', 'vref_grid', 'ccrs']
diff --git a/tests/test_demcollection.py b/tests/test_demcollection.py
index b14df949..1f432439 100644
--- a/tests/test_demcollection.py
+++ b/tests/test_demcollection.py
@@ -24,10 +24,10 @@ def test_init(self) -> None:
         # Make sure the glacier was bigger in 1990, since this is assumed later.
         assert scott_1990.ds.area.sum() > scott_2010.ds.area.sum()
 
-        mask_2010 = scott_2010.create_mask(self.dem_2009).reshape(self.dem_2009.data.shape)
+        mask_2010 = scott_2010.create_mask(self.dem_2009)
 
         dem_2060 = self.dem_2009.copy()
-        dem_2060.data[mask_2010] -= 30
+        dem_2060[mask_2010] -= 30
 
         dems = xdem.DEMCollection(
             [self.dem_1990, self.dem_2009, dem_2060],
@@ -36,7 +36,7 @@ def test_init(self) -> None:
             reference_dem=1,
         )
 
-        # Check that the first raster is the oldest one and
+        # Check that the first raster is the oldest one
         assert dems.dems[0].data.max() == self.dem_1990.data.max()
         assert dems.reference_dem.data.max() == self.dem_2009.data.max()
 
diff --git a/tests/test_examples.py b/tests/test_examples.py
index f6af544b..f6e028b4 100644
--- a/tests/test_examples.py
+++ b/tests/test_examples.py
@@ -62,6 +62,6 @@ def test_array_nodata(self, rst_and_truenodata: tuple[Raster, int]) -> None:
 
         rst = rst_and_truenodata[0]
         truenodata = rst_and_truenodata[1]
-        mask = gu.spatial_tools.get_array_and_mask(rst)[1]
+        mask = gu.raster.get_array_and_mask(rst)[1]
 
         assert np.sum(mask) == truenodata
diff --git a/tests/test_filters.py b/tests/test_filters.py
index 3cc00e33..30c3c659 100644
--- a/tests/test_filters.py
+++ b/tests/test_filters.py
@@ -12,8 +12,8 @@ class TestFilters:
     """Test cases for the filter functions."""
 
     # Load example data.
-    dem_2009 = gu.georaster.Raster(xdem.examples.get_path("longyearbyen_ref_dem"))
-    dem_1990 = gu.georaster.Raster(xdem.examples.get_path("longyearbyen_tba_dem")).reproject(dem_2009, silent=True)
+    dem_2009 = gu.Raster(xdem.examples.get_path("longyearbyen_ref_dem"))
+    dem_1990 = gu.Raster(xdem.examples.get_path("longyearbyen_tba_dem")).reproject(dem_2009, silent=True)
 
     def test_gauss(self) -> None:
         """Test applying the various Gaussian filters on DEMs with/without NaNs"""
@@ -51,7 +51,7 @@ def test_gauss(self) -> None:
         assert np.nanmax(dem_with_nans) > np.max(dem_sm)
 
         # Test that it works with 3D arrays
-        array_3d = np.vstack((dem_array, dem_array))
+        array_3d = np.vstack((dem_array[np.newaxis, :], dem_array[np.newaxis, :]))
         dem_sm = xdem.filters.gaussian_filter_scipy(array_3d, sigma=5)
         assert array_3d.shape == dem_sm.shape
 
@@ -59,7 +59,7 @@ def test_gauss(self) -> None:
         pytest.raises(NotImplementedError, xdem.filters.gaussian_filter_cv, array_3d, sigma=5)
 
         # Tests that it fails with 1D arrays with appropriate error
-        data = dem_array[0, :, 0]
+        data = dem_array[:, 0]
         pytest.raises(ValueError, xdem.filters.gaussian_filter_scipy, data, sigma=5)
         pytest.raises(ValueError, xdem.filters.gaussian_filter_cv, data, sigma=5)
 
@@ -73,19 +73,19 @@ def test_dist_filter(self) -> None:
         count = 1000
         cols = np.random.randint(0, high=self.dem_1990.width - 1, size=count, dtype=int)
         rows = np.random.randint(0, high=self.dem_1990.height - 1, size=count, dtype=int)
-        ddem.data[0, rows, cols] = 5000
+        ddem.data[rows, cols] = 5000
 
         # Filter gross outliers
         filtered_ddem = xdem.filters.distance_filter(ddem.data, radius=20, outlier_threshold=50)
 
         # Check that all outliers were properly filtered
-        assert np.all(np.isnan(filtered_ddem[0, rows, cols]))
+        assert np.all(np.isnan(filtered_ddem[rows, cols]))
 
         # Assert that non filtered pixels remain the same
         assert ddem.data.shape == filtered_ddem.shape
         assert np.all(ddem.data[np.isfinite(filtered_ddem)] == filtered_ddem[np.isfinite(filtered_ddem)])
 
         # Check that it works with NaNs too
-        ddem.data[0, rows[:500], cols[:500]] = np.nan
+        ddem.data[rows[:500], cols[:500]] = np.nan
         filtered_ddem = xdem.filters.distance_filter(ddem.data, radius=20, outlier_threshold=50)
-        assert np.all(np.isnan(filtered_ddem[0, rows, cols]))
+        assert np.all(np.isnan(filtered_ddem[rows, cols]))
diff --git a/tests/test_misc.py b/tests/test_misc.py
index 4ba10102..ab76b073 100644
--- a/tests/test_misc.py
+++ b/tests/test_misc.py
@@ -105,3 +105,54 @@ def useless_func() -> int:
         else:
             with pytest.raises(ValueError, match="^" + text + "$"):
                 useless_func()
+
+    def test_diff_environment_yml(self, capsys) -> None:  # type: ignore
+
+        # Test with synthetic environment
+        env = {"dependencies": ["python==3.9", "numpy", "fiona"]}
+        devenv = {"dependencies": ["python==3.9", "numpy", "fiona", "opencv"]}
+
+        # This should print the difference between the two
+        xdem.misc.diff_environment_yml(env, devenv, input_dict=True, print_dep="conda")
+
+        # Capture the stdout and check it is indeed the right diff
+        captured = capsys.readouterr().out
+        assert captured == "opencv\n"
+
+        # This should print the difference including pip
+        xdem.misc.diff_environment_yml(env, devenv, input_dict=True, print_dep="both")
+
+        captured = capsys.readouterr().out
+        assert captured == "opencv\nNone\n"
+
+        env2 = {"dependencies": ["python==3.9", "numpy", "fiona"]}
+        devenv2 = {"dependencies": ["python==3.9", "numpy", "fiona", "opencv", {"pip": ["geoutils", "-e ./"]}]}
+
+        # The diff function should not account for -e ./ that is the local install for developers
+        xdem.misc.diff_environment_yml(env2, devenv2, input_dict=True, print_dep="both")
+        captured = capsys.readouterr().out
+
+        assert captured == "opencv\ngeoutils\n"
+
+        # This should print only pip
+        xdem.misc.diff_environment_yml(env2, devenv2, input_dict=True, print_dep="pip")
+        captured = capsys.readouterr().out
+
+        assert captured == "geoutils\n"
+
+        # This should raise an error because print_dep is not well defined
+        with pytest.raises(ValueError, match='The argument "print_dep" can only be "conda", "pip" or "both".'):
+            xdem.misc.diff_environment_yml(env2, devenv2, input_dict=True, print_dep="lol")
+
+        # When the dependencies are not defined in dev-env but in env, it should raise an error
+        # For normal dependencies
+        env3 = {"dependencies": ["python==3.9", "numpy", "fiona", "lol"]}
+        devenv3 = {"dependencies": ["python==3.9", "numpy", "fiona", "opencv", {"pip": ["geoutils"]}]}
+        with pytest.raises(ValueError, match="The following dependencies are listed in env but not dev-env: lol"):
+            xdem.misc.diff_environment_yml(env3, devenv3, input_dict=True, print_dep="pip")
+
+        # For pip dependencies
+        env4 = {"dependencies": ["python==3.9", "numpy", "fiona", {"pip": ["lol"]}]}
+        devenv4 = {"dependencies": ["python==3.9", "numpy", "fiona", "opencv", {"pip": ["geoutils"]}]}
+        with pytest.raises(ValueError, match="The following pip dependencies are listed in env but not dev-env: lol"):
+            xdem.misc.diff_environment_yml(env4, devenv4, input_dict=True, print_dep="pip")
diff --git a/tests/test_spatialstats.py b/tests/test_spatialstats.py
index 07107b49..3f89915e 100644
--- a/tests/test_spatialstats.py
+++ b/tests/test_spatialstats.py
@@ -16,7 +16,7 @@
 import xdem
 from xdem import examples
 from xdem._typing import NDArrayf
-from xdem.spatialstats import EmpiricalVariogramKArgs
+from xdem.spatialstats import EmpiricalVariogramKArgs, nmad
 
 with warnings.catch_warnings():
     warnings.filterwarnings("ignore", category=DeprecationWarning)
@@ -37,6 +37,28 @@ def load_ref_and_diff() -> tuple[Raster, Raster, NDArrayf, Vector]:
     return reference_raster, ddem, mask, outlines
 
 
+class TestStats:
+
+    # Load data for the entire test class
+    ref, diff, mask, outlines = load_ref_and_diff()
+
+    def test_nmad(self) -> None:
+        """Test NMAD functionality runs on any type of input"""
+
+        # Check that the NMAD is computed the same with a raster, masked array or NaN array
+        nmad_raster = nmad(self.diff)
+        nmad_ma = nmad(self.diff.data)
+        nmad_array = nmad(self.diff.get_nanarray())
+
+        assert nmad_raster == nmad_ma == nmad_array
+
+        # Check that the scaling factor works
+        nmad_1 = nmad(self.diff, nfact=1)
+        nmad_2 = nmad(self.diff, nfact=2)
+
+        assert nmad_1 * 2 == nmad_2
+
+
 class TestBinning:
 
     # Load data for the entire test class
@@ -51,7 +73,7 @@ def test_nd_binning(self) -> None:
         """Check that the nd_binning function works adequately and save dataframes to files for later tests"""
 
         # Subsampler
-        indices = gu.spatial_tools.subsample_raster(
+        indices = gu.raster.subsample_array(
             self.diff.data.flatten(), subsample=10000, return_indices=True, random_state=42
         )
 
@@ -259,16 +281,10 @@ def test_interp_nd_binning_realdata(self) -> None:
     def test_two_step_standardization(self) -> None:
         """Test two-step standardization function"""
 
-        # Test this gives the same results as when using the base functions
-        diff_arr = gu.spatial_tools.get_array_and_mask(self.diff)[0]
-        slope_arr = gu.spatial_tools.get_array_and_mask(self.slope)[0]
-        maximum_curv_arr = gu.spatial_tools.get_array_and_mask(self.maximum_curv)[0]
-        stable_mask_arr = ~self.outlines.create_mask(self.ref).squeeze()
-
         # Reproduce the first steps of binning
         df_binning = xdem.spatialstats.nd_binning(
-            values=diff_arr[stable_mask_arr],
-            list_var=[slope_arr[stable_mask_arr], maximum_curv_arr[stable_mask_arr]],
+            values=self.diff[~self.mask],
+            list_var=[self.slope[~self.mask], self.maximum_curv[~self.mask]],
             list_var_names=["var1", "var2"],
             statistics=[xdem.spatialstats.nmad],
         )
@@ -276,9 +292,7 @@ def test_two_step_standardization(self) -> None:
             df_binning, list_var_names=["var1", "var2"], statistic="nmad"
         )
         # The zscore spread should not be one right after binning
-        zscores = diff_arr[stable_mask_arr] / unscaled_fun(
-            (slope_arr[stable_mask_arr], maximum_curv_arr[stable_mask_arr])
-        )
+        zscores = self.diff[~self.mask] / unscaled_fun((self.slope[~self.mask], self.maximum_curv[~self.mask]))
         scale_fac = xdem.spatialstats.nmad(zscores)
         assert scale_fac != 1
 
@@ -288,8 +302,8 @@ def test_two_step_standardization(self) -> None:
         scale_fac_std = np.nanstd(zscores)
         zscores /= scale_fac_std
         zscores_2, final_func = xdem.spatialstats.two_step_standardization(
-            diff_arr[stable_mask_arr],
-            list_var=[slope_arr[stable_mask_arr], maximum_curv_arr[stable_mask_arr]],
+            dvalues=self.diff[~self.mask],
+            list_var=[self.slope[~self.mask], self.maximum_curv[~self.mask]],
             unscaled_error_fun=unscaled_fun,
             spread_statistic=np.nanstd,
             fac_spread_outliers=3,
@@ -311,15 +325,9 @@ def test_estimate_model_heteroscedasticity_and_infer_from_stable(self) -> None:
             dvalues=self.diff, list_var=[self.slope, self.maximum_curv], unstable_mask=self.outlines
         )
 
-        # Test this gives the same results as when using the base functions
-        diff_arr = gu.spatial_tools.get_array_and_mask(self.diff)[0]
-        slope_arr = gu.spatial_tools.get_array_and_mask(self.slope)[0]
-        maximum_curv_arr = gu.spatial_tools.get_array_and_mask(self.maximum_curv)[0]
-        stable_mask_arr = ~self.outlines.create_mask(self.ref).squeeze()
-
         df_binning_2, err_fun_2 = xdem.spatialstats.estimate_model_heteroscedasticity(
-            dvalues=diff_arr[stable_mask_arr],
-            list_var=[slope_arr[stable_mask_arr], maximum_curv_arr[stable_mask_arr]],
+            dvalues=self.diff[~self.mask],
+            list_var=[self.slope[~self.mask], self.maximum_curv[~self.mask]],
             list_var_names=["var1", "var2"],
         )
 
@@ -329,8 +337,8 @@ def test_estimate_model_heteroscedasticity_and_infer_from_stable(self) -> None:
         assert np.array_equal(err_fun_1((test_slopes, test_max_curvs)), err_fun_2((test_slopes, test_max_curvs)))
 
         # Test the error map is consistent as well
-        errors_2_arr = err_fun_2((slope_arr, maximum_curv_arr))
-        errors_1_arr = gu.spatial_tools.get_array_and_mask(errors_1)[0]
+        errors_2_arr = err_fun_2((self.slope.get_nanarray(), self.maximum_curv.get_nanarray()))
+        errors_1_arr = gu.raster.get_array_and_mask(errors_1)[0]
         assert np.array_equal(errors_1_arr, errors_2_arr, equal_nan=True)
 
         # Save for use in TestVariogram
@@ -339,15 +347,15 @@ def test_estimate_model_heteroscedasticity_and_infer_from_stable(self) -> None:
         # Check that errors are raised with wrong input
         with pytest.raises(ValueError, match="The values must be a Raster or NumPy array, or a list of those."):
             xdem.spatialstats.infer_heteroscedasticity_from_stable(
-                dvalues="not_an_array", stable_mask=~self.mask.squeeze(), list_var=[slope_arr]
+                dvalues="not_an_array", stable_mask=~self.mask, list_var=[self.slope.get_nanarray()]
             )
-        with pytest.raises(ValueError, match="The stable mask must be a Vector, GeoDataFrame or NumPy array."):
+        with pytest.raises(ValueError, match="The stable mask must be a Vector, Mask, GeoDataFrame or NumPy array."):
             xdem.spatialstats.infer_heteroscedasticity_from_stable(
-                dvalues=self.diff, stable_mask="not_a_vector_or_array", list_var=[slope_arr]
+                dvalues=self.diff, stable_mask="not_a_vector_or_array", list_var=[self.slope.get_nanarray()]
             )
-        with pytest.raises(ValueError, match="The unstable mask must be a Vector, GeoDataFrame or NumPy array."):
+        with pytest.raises(ValueError, match="The unstable mask must be a Vector, Mask, GeoDataFrame or NumPy array."):
             xdem.spatialstats.infer_heteroscedasticity_from_stable(
-                dvalues=self.diff, unstable_mask="not_a_vector_or_array", list_var=[slope_arr]
+                dvalues=self.diff, unstable_mask="not_a_vector_or_array", list_var=[self.slope.get_nanarray()]
             )
 
         with pytest.raises(
@@ -356,7 +364,7 @@ def test_estimate_model_heteroscedasticity_and_infer_from_stable(self) -> None:
             "values contain a Raster.",
         ):
             xdem.spatialstats.infer_heteroscedasticity_from_stable(
-                dvalues=diff_arr, stable_mask=self.outlines, list_var=[slope_arr]
+                dvalues=self.diff.get_nanarray(), stable_mask=self.outlines, list_var=[self.slope.get_nanarray()]
             )
 
     def test_plot_binning(self) -> None:
@@ -466,7 +474,7 @@ def test_sample_empirical_variogram_speed(self) -> None:
         )
 
         # Index of valid values
-        values_arr, mask_nodata = gu.spatial_tools.get_array_and_mask(values)
+        values_arr, mask_nodata = gu.raster.get_array_and_mask(values)
 
         t3 = time.time()
         rems = skgstat.RasterEquidistantMetricSpace(
@@ -759,17 +767,17 @@ def test_estimate_model_spatial_correlation_and_infer_from_stable(self) -> None:
 
         # Run wrapper infer from stable function with a Raster and the mask, and check the consistency there as well
         emp_vgm_3, params_model_vgm_3, _ = xdem.spatialstats.infer_spatial_correlation_from_stable(
-            dvalues=zscores, stable_mask=~self.mask.squeeze(), list_models=["Gau", "Sph"], subsample=10, random_state=42
+            dvalues=zscores, stable_mask=~self.mask, list_models=["Gau", "Sph"], subsample=10, random_state=42
         )
         pd.testing.assert_frame_equal(emp_vgm_1, emp_vgm_3)
         pd.testing.assert_frame_equal(params_model_vgm_1, params_model_vgm_3)
 
         # Run again with array instead of Raster as input
-        zscores_arr = gu.spatial_tools.get_array_and_mask(zscores)[0]
+        zscores_arr = gu.raster.get_array_and_mask(zscores)[0]
         emp_vgm_4, params_model_vgm_4, _ = xdem.spatialstats.infer_spatial_correlation_from_stable(
             dvalues=zscores_arr,
             gsd=self.diff.res[0],
-            stable_mask=~self.mask.squeeze(),
+            stable_mask=~self.mask,
             list_models=["Gau", "Sph"],
             subsample=10,
             random_state=42,
@@ -781,7 +789,7 @@ def test_estimate_model_spatial_correlation_and_infer_from_stable(self) -> None:
         _, params_model_vgm_5, _ = xdem.spatialstats.infer_spatial_correlation_from_stable(
             dvalues=zscores_arr,
             gsd=self.diff.res[0],
-            stable_mask=~self.mask.squeeze(),
+            stable_mask=~self.mask,
             list_models=["Gau", "Sph"],
             subsample=200,
             random_state=42,
@@ -794,17 +802,17 @@ def test_estimate_model_spatial_correlation_and_infer_from_stable(self) -> None:
         # Check that errors are raised with wrong input
         with pytest.raises(ValueError, match="The values must be a Raster or NumPy array, or a list of those."):
             xdem.spatialstats.infer_spatial_correlation_from_stable(
-                dvalues="not_an_array", stable_mask=~self.mask.squeeze(), list_models=["Gau", "Sph"], random_state=42
+                dvalues="not_an_array", stable_mask=~self.mask, list_models=["Gau", "Sph"], random_state=42
             )
-        with pytest.raises(ValueError, match="The stable mask must be a Vector, GeoDataFrame or NumPy array."):
+        with pytest.raises(ValueError, match="The stable mask must be a Vector, Mask, GeoDataFrame or NumPy array."):
             xdem.spatialstats.infer_spatial_correlation_from_stable(
                 dvalues=self.diff, stable_mask="not_a_vector_or_array", list_models=["Gau", "Sph"], random_state=42
             )
-        with pytest.raises(ValueError, match="The unstable mask must be a Vector, GeoDataFrame or NumPy array."):
+        with pytest.raises(ValueError, match="The unstable mask must be a Vector, Mask, GeoDataFrame or NumPy array."):
             xdem.spatialstats.infer_spatial_correlation_from_stable(
                 dvalues=self.diff, unstable_mask="not_a_vector_or_array", list_models=["Gau", "Sph"], random_state=42
             )
-        diff_on_stable_arr = gu.spatial_tools.get_array_and_mask(diff_on_stable)[0]
+        diff_on_stable_arr = gu.raster.get_array_and_mask(diff_on_stable)[0]
         with pytest.raises(
             ValueError,
             match="The stable mask can only passed as a Vector or GeoDataFrame if the input "
@@ -1008,7 +1016,7 @@ def test_number_effective_samples(self) -> None:
             subsample=10,
         )
         # Second, get coordinates manually and compute with the neff_approx_hugonnet function
-        mask = outlines_brom.create_mask(xres=res)
+        mask = outlines_brom.create_mask(xres=res, as_array=True)
         x = res * np.arange(0, mask.shape[0])
         y = res * np.arange(0, mask.shape[1])
         coords = np.array(np.meshgrid(y, x))
@@ -1166,7 +1174,7 @@ def test_patches_method_loop_quadrant(self) -> None:
         # Check the patches method runs
         df, df_full = xdem.spatialstats.patches_method(
             diff,
-            unstable_mask=mask.squeeze(),
+            unstable_mask=mask,
             gsd=gsd,
             areas=[area],
             random_state=42,
@@ -1205,7 +1213,7 @@ def test_patches_method_convolution(self) -> None:
         # First, the patches method runs with scipy
         df = xdem.spatialstats.patches_method(
             diff,
-            unstable_mask=mask.squeeze(),
+            unstable_mask=mask,
             gsd=gsd,
             areas=[area, area * 10],
             random_state=42,
diff --git a/tests/test_terrain.py b/tests/test_terrain.py
index 124fa221..2c7e7a99 100644
--- a/tests/test_terrain.py
+++ b/tests/test_terrain.py
@@ -208,8 +208,8 @@ def test_attribute_functions_against_richdem(self, attribute: str) -> None:
         dem = self.dem.copy()
 
         # Derive the attribute using both RichDEM and xdem
-        attr_xdem = gu.spatial_tools.get_array_and_mask(functions_xdem[attribute](dem))[0].squeeze()
-        attr_richdem = gu.spatial_tools.get_array_and_mask(functions_richdem[attribute](dem))[0].squeeze()
+        attr_xdem = gu.raster.get_array_and_mask(functions_xdem[attribute](dem))[0].squeeze()
+        attr_richdem = gu.raster.get_array_and_mask(functions_richdem[attribute](dem))[0].squeeze()
 
         # We compute the difference and keep only valid values
         diff = attr_xdem - attr_richdem
diff --git a/tests/test_volume.py b/tests/test_volume.py
index 032c4dfa..338f8198 100644
--- a/tests/test_volume.py
+++ b/tests/test_volume.py
@@ -13,9 +13,9 @@ class TestLocalHypsometric:
     """Test cases for the local hypsometric method."""
 
     # Load example data.
-    dem_2009 = gu.georaster.Raster(xdem.examples.get_path("longyearbyen_ref_dem"))
-    dem_1990 = gu.georaster.Raster(xdem.examples.get_path("longyearbyen_tba_dem")).reproject(dem_2009, silent=True)
-    outlines = gu.geovector.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines"))
+    dem_2009 = gu.Raster(xdem.examples.get_path("longyearbyen_ref_dem"))
+    dem_1990 = gu.Raster(xdem.examples.get_path("longyearbyen_tba_dem")).reproject(dem_2009, silent=True)
+    outlines = gu.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines"))
     all_outlines = outlines.copy()
 
     # Filter to only look at the Scott Turnerbreen glacier
@@ -25,27 +25,26 @@ class TestLocalHypsometric:
 
     def test_bin_ddem(self) -> None:
         """Test dDEM binning."""
-        ddem = self.dem_2009.data - self.dem_1990.data
+        ddem = self.dem_2009 - self.dem_1990
 
-        ddem_bins = xdem.volume.hypsometric_binning(
-            ddem[self.mask], self.dem_2009.data[self.mask], bins=50, kind="fixed"
-        )
+        ddem_bins = xdem.volume.hypsometric_binning(ddem[self.mask], self.dem_2009[self.mask], bins=50, kind="fixed")
 
         ddem_bins_masked = xdem.volume.hypsometric_binning(
-            np.ma.masked_array(ddem, mask=~self.mask), np.ma.masked_array(self.dem_2009.data, mask=~self.mask)
+            np.ma.masked_array(ddem.data, mask=~self.mask.data.filled(False)),
+            np.ma.masked_array(self.dem_2009.data, mask=~self.mask.data.filled(False)),
         )
 
         ddem_stds = xdem.volume.hypsometric_binning(
-            ddem[self.mask], self.dem_2009.data[self.mask], aggregation_function=np.std
+            ddem[self.mask], self.dem_2009[self.mask], aggregation_function=np.std
         )
         assert ddem_stds["value"].mean() < 50
         assert np.abs(np.mean(ddem_bins["value"] - ddem_bins_masked["value"])) < 0.01
 
     def test_interpolate_ddem_bins(self) -> pd.Series:
         """Test dDEM bin interpolation."""
-        ddem = self.dem_2009.data - self.dem_1990.data
+        ddem = self.dem_2009 - self.dem_1990
 
-        ddem_bins = xdem.volume.hypsometric_binning(ddem[self.mask], self.dem_2009.data[self.mask])
+        ddem_bins = xdem.volume.hypsometric_binning(ddem[self.mask], self.dem_2009[self.mask])
 
         # Simulate a missing bin
         ddem_bins.iloc[3, 0] = np.nan
@@ -71,17 +70,17 @@ def test_area_calculation(self) -> None:
 
         # Test the area calculation with normal parameters.
         bin_area = xdem.volume.calculate_hypsometry_area(
-            ddem_bins, self.dem_2009.data[self.mask], pixel_size=self.dem_2009.res[0]
+            ddem_bins, self.dem_2009[self.mask], pixel_size=self.dem_2009.res[0]
         )
 
         # Test area calculation with differing pixel x/y resolution.
         xdem.volume.calculate_hypsometry_area(
-            ddem_bins, self.dem_2009.data[self.mask], pixel_size=(self.dem_2009.res[0], self.dem_2009.res[0] + 1)
+            ddem_bins, self.dem_2009[self.mask], pixel_size=(self.dem_2009.res[0], self.dem_2009.res[0] + 1)
         )
 
         # Add some nans to the reference DEM
-        data_with_nans = self.dem_2009.data[self.mask]
-        data_with_nans[[2, 5]] = np.nan
+        data_with_nans = self.dem_2009
+        data_with_nans.data[2, 5] = np.nan
 
         # Make sure that the above results in the correct error.
         try:
@@ -93,7 +92,7 @@ def test_area_calculation(self) -> None:
         # Try to pass an incorrect timeframe= parameter
         try:
             xdem.volume.calculate_hypsometry_area(
-                ddem_bins, self.dem_2009.data[self.mask], pixel_size=self.dem_2009.res[0], timeframe="blergh"
+                ddem_bins, self.dem_2009[self.mask], pixel_size=self.dem_2009.res[0], timeframe="blergh"
             )
         except ValueError as exception:
             if "Argument 'timeframe=blergh' is invalid" not in str(exception):
@@ -102,9 +101,7 @@ def test_area_calculation(self) -> None:
         # Mess up the dDEM bins and see if it gives the correct error
         ddem_bins.iloc[3] = np.nan
         try:
-            xdem.volume.calculate_hypsometry_area(
-                ddem_bins, self.dem_2009.data[self.mask], pixel_size=self.dem_2009.res[0]
-            )
+            xdem.volume.calculate_hypsometry_area(ddem_bins, self.dem_2009[self.mask], pixel_size=self.dem_2009.res[0])
         except AssertionError as exception:
             if "cannot contain NaNs" not in str(exception):
                 raise exception
@@ -114,19 +111,19 @@ def test_area_calculation(self) -> None:
 
     def test_ddem_bin_methods(self) -> None:
         """Test different dDEM binning methods."""
-        ddem = self.dem_2009.data - self.dem_1990.data
+        ddem = self.dem_2009 - self.dem_1990
 
         # equal height is already tested in test_bin_ddem
 
         # Make a fixed amount of bins
         equal_count_bins = xdem.volume.hypsometric_binning(
-            ddem[self.mask], self.dem_2009.data[self.mask], bins=50, kind="count"
+            ddem[self.mask], self.dem_2009[self.mask], bins=50, kind="count"
         )
         assert equal_count_bins.shape[0] == 50
 
         # Make 50 bins with approximately the same area (pixel count)
         quantile_bins = xdem.volume.hypsometric_binning(
-            ddem[self.mask], self.dem_2009.data[self.mask], bins=50, kind="quantile"
+            ddem[self.mask], self.dem_2009[self.mask], bins=50, kind="quantile"
         )
 
         assert quantile_bins.shape[0] == 50
@@ -136,7 +133,7 @@ def test_ddem_bin_methods(self) -> None:
         # Try to feed the previous bins as "custom" bins to the function.
         custom_bins = xdem.volume.hypsometric_binning(
             ddem[self.mask],
-            self.dem_2009.data[self.mask],
+            self.dem_2009[self.mask],
             bins=np.r_[quantile_bins.index.left[0], quantile_bins.index.right],
             kind="custom",
         )
@@ -145,9 +142,9 @@ def test_ddem_bin_methods(self) -> None:
 
 
 class TestNormHypsometric:
-    dem_2009 = gu.georaster.Raster(xdem.examples.get_path("longyearbyen_ref_dem"))
-    dem_1990 = gu.georaster.Raster(xdem.examples.get_path("longyearbyen_tba_dem")).reproject(dem_2009, silent=True)
-    outlines = gu.geovector.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines"))
+    dem_2009 = gu.Raster(xdem.examples.get_path("longyearbyen_ref_dem"))
+    dem_1990 = gu.Raster(xdem.examples.get_path("longyearbyen_tba_dem")).reproject(dem_2009, silent=True)
+    outlines = gu.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines"))
 
     glacier_index_map = outlines.rasterize(dem_2009)
     ddem = dem_2009.data - dem_1990.data
@@ -157,7 +154,7 @@ def test_regional_signal(self, n_bins: int) -> None:
         warnings.simplefilter("error")
 
         signal = xdem.volume.get_regional_hypsometric_signal(
-            ddem=self.ddem, ref_dem=self.dem_2009.data, glacier_index_map=self.glacier_index_map, n_bins=n_bins
+            ddem=self.ddem, ref_dem=self.dem_2009, glacier_index_map=self.glacier_index_map, n_bins=n_bins
         )
 
         assert signal["w_mean"].min() >= 0.0
@@ -212,10 +209,10 @@ def test_regional_hypsometric_interp(self) -> None:
         warnings.simplefilter("error")
 
         # Extract a normalized regional hypsometric signal.
-        ddem = self.dem_2009.data - self.dem_1990.data
+        ddem = self.dem_2009 - self.dem_1990
 
         signal = xdem.volume.get_regional_hypsometric_signal(
-            ddem=self.ddem, ref_dem=self.dem_2009.data, glacier_index_map=self.glacier_index_map
+            ddem=self.ddem, ref_dem=self.dem_2009, glacier_index_map=self.glacier_index_map
         )
 
         if False:
@@ -235,19 +232,19 @@ def test_regional_hypsometric_interp(self) -> None:
 
         # Try the normalized regional hypsometric interpolation.
         # Synthesize random nans in 80% of the data.
-        ddem.mask.ravel()[np.random.choice(ddem.data.size, int(ddem.data.size * 0.80), replace=False)] = True
+        ddem.data.mask.ravel()[np.random.choice(ddem.data.size, int(ddem.data.size * 0.80), replace=False)] = True
         # Fill the dDEM using the de-normalized signal.
         filled_ddem = xdem.volume.norm_regional_hypsometric_interpolation(
-            voided_ddem=ddem, ref_dem=self.dem_2009.data, glacier_index_map=self.glacier_index_map
+            voided_ddem=ddem, ref_dem=self.dem_2009, glacier_index_map=self.glacier_index_map
         )
         # Fill the dDEM using the de-normalized signal and create an idealized dDEM
         idealized_ddem = xdem.volume.norm_regional_hypsometric_interpolation(
-            voided_ddem=ddem, ref_dem=self.dem_2009.data, glacier_index_map=self.glacier_index_map, idealized_ddem=True
+            voided_ddem=ddem, ref_dem=self.dem_2009, glacier_index_map=self.glacier_index_map, idealized_ddem=True
         )
         assert not np.array_equal(filled_ddem, idealized_ddem)
 
         # Check that all glacier-values are finite
-        assert np.count_nonzero(np.isnan(idealized_ddem)[self.glacier_index_map > 0]) == 0
+        assert np.count_nonzero(np.isnan(idealized_ddem)[self.glacier_index_map.get_nanarray() > 0]) == 0
         # Validate that the un-idealized dDEM has a higher gradient variance (more ups and downs)
         filled_gradient = np.linalg.norm(np.gradient(filled_ddem), axis=0)
         ideal_gradient = np.linalg.norm(np.gradient(idealized_ddem), axis=0)
@@ -264,9 +261,9 @@ def test_regional_hypsometric_interp(self) -> None:
             plt.show()
 
         # Extract the finite glacier values.
-        changes = ddem.data.squeeze()[self.glacier_index_map > 0]
+        changes = ddem.data.squeeze()[self.glacier_index_map.get_nanarray() > 0]
         changes = changes[np.isfinite(changes)]
-        interp_changes = filled_ddem[self.glacier_index_map > 0]
+        interp_changes = filled_ddem[self.glacier_index_map.get_nanarray() > 0]
         interp_changes = interp_changes[np.isfinite(interp_changes)]
 
         # Validate that the interpolated (20% data) means and stds are similar to the original (100% data)
diff --git a/xdem/coreg.py b/xdem/coreg.py
index 6694cf09..a032e5c1 100644
--- a/xdem/coreg.py
+++ b/xdem/coreg.py
@@ -25,9 +25,15 @@
 import scipy.ndimage
 import scipy.optimize
 import skimage.transform
-from geoutils import spatial_tools
 from geoutils._typing import AnyNumber
-from geoutils.georaster import RasterType, raster
+from geoutils.raster import (
+    Mask,
+    RasterType,
+    get_array_and_mask,
+    raster,
+    subdivide_array,
+    subsample_array,
+)
 from rasterio import Affine
 from tqdm import tqdm, trange
 
@@ -293,7 +299,7 @@ def residuals(coefs: NDArrayf, x_coords: NDArrayf, y_coords: NDArrayf, targets:
         print("Estimating deramp function...")
 
     # reduce number of elements for speed
-    rand_indices = gu.spatial_tools.subsample_raster(x_coords, subsample=subsample, return_indices=True)
+    rand_indices = subsample_array(x_coords, subsample=subsample, return_indices=True)
     x_coords = x_coords[rand_indices]
     y_coords = y_coords[rand_indices]
     ddem = ddem[rand_indices]
@@ -319,9 +325,7 @@ def fit_ramp(x: NDArrayf, y: NDArrayf) -> NDArrayf:
     return fit_ramp, coefs
 
 
-def mask_as_array(
-    reference_raster: gu.georaster.Raster, mask: str | gu.geovector.Vector | gu.georaster.Raster
-) -> NDArrayf:
+def mask_as_array(reference_raster: gu.Raster, mask: str | gu.Vector | gu.Raster) -> NDArrayf:
     """
     Convert a given mask into an array.
 
@@ -337,27 +341,26 @@ def mask_as_array(
     if isinstance(mask, str):
         # First try to load it as a Vector
         try:
-            mask = gu.geovector.Vector(mask)
+            mask = gu.Vector(mask)
         # If the format is unsopported, try loading as a Raster
         except fiona.errors.DriverError:
             try:
-                mask = gu.georaster.Raster(mask)
+                mask = gu.Raster(mask)
             # If that fails, raise an error
             except rio.errors.RasterioIOError:
                 raise ValueError(f"Mask path not in a supported Raster or Vector format: {mask}")
 
     # At this point, the mask variable is either a Raster or a Vector
     # Now, convert the mask into an array by either rasterizing a Vector or by fetching a Raster's data
-    if isinstance(mask, gu.geovector.Vector):
-        mask_array = mask.create_mask(reference_raster)
-    elif isinstance(mask, gu.georaster.Raster):
+    if isinstance(mask, gu.Vector):
+        mask_array = mask.create_mask(reference_raster, as_array=True)
+    elif isinstance(mask, gu.Raster):
         # The true value is the maximum value in the raster, unless the maximum value is 0 or False
         true_value = np.nanmax(mask.data) if not np.nanmax(mask.data) in [0, False] else True
         mask_array = (mask.data == true_value).squeeze()
     else:
         raise TypeError(
-            f"Mask has invalid type: {type(mask)}. Expected one of: "
-            f"{[gu.georaster.Raster, gu.geovector.Vector, str, type(None)]}"
+            f"Mask has invalid type: {type(mask)}. Expected one of: " f"{[gu.Raster, gu.Vector, str, type(None)]}"
         )
 
     return mask_array
@@ -433,7 +436,7 @@ def fit(
         self: CoregType,
         reference_dem: NDArrayf | MArrayf | RasterType,
         dem_to_be_aligned: NDArrayf | MArrayf | RasterType,
-        inlier_mask: NDArrayf | None = None,
+        inlier_mask: NDArrayf | Mask | None = None,
         transform: rio.transform.Affine | None = None,
         crs: rio.crs.CRS | None = None,
         weights: NDArrayf | None = None,
@@ -496,13 +499,16 @@ def fit(
         if crs is None:
             raise ValueError("'crs' must be given if both DEMs are array-like.")
 
-        ref_dem, ref_mask = spatial_tools.get_array_and_mask(reference_dem)
-        tba_dem, tba_mask = spatial_tools.get_array_and_mask(dem_to_be_aligned)
+        ref_dem, ref_mask = get_array_and_mask(reference_dem)
+        tba_dem, tba_mask = get_array_and_mask(dem_to_be_aligned)
 
         # Make sure that the mask has an expected format.
         if inlier_mask is not None:
-            inlier_mask = np.asarray(inlier_mask).squeeze()
-            assert inlier_mask.dtype == bool, f"Invalid mask dtype: '{inlier_mask.dtype}'. Expected 'bool'"
+            if isinstance(inlier_mask, Mask):
+                inlier_mask = inlier_mask.data.filled(False).squeeze()
+            else:
+                inlier_mask = np.asarray(inlier_mask).squeeze()
+                assert inlier_mask.dtype == bool, f"Invalid mask dtype: '{inlier_mask.dtype}'. Expected 'bool'"
 
             if np.all(~inlier_mask):
                 raise ValueError("'inlier_mask' had no inliers.")
@@ -521,7 +527,7 @@ def fit(
             full_mask = (
                 ~ref_mask & ~tba_mask & (np.asarray(inlier_mask) if inlier_mask is not None else True)
             ).squeeze()
-            random_indices = gu.spatial_tools.subsample_raster(full_mask, subsample=subsample, return_indices=True)
+            random_indices = subsample_array(full_mask, subsample=subsample, return_indices=True)
             full_mask[random_indices] = False
 
         # Run the associated fitting function
@@ -607,7 +613,7 @@ def apply(
                 raise ValueError("'crs' must be given if DEM is array-like.")
 
         # The array to provide the functions will be an ndarray with NaNs for masked out areas.
-        dem_array, dem_mask = spatial_tools.get_array_and_mask(dem)
+        dem_array, dem_mask = get_array_and_mask(dem)
 
         if np.all(dem_mask):
             raise ValueError("'dem' had only NaNs")
@@ -776,7 +782,7 @@ def residuals(
         aligned_dem, _ = self.apply(dem_to_be_aligned, transform=transform, crs=crs)
 
         # Format the reference DEM
-        ref_arr, ref_mask = spatial_tools.get_array_and_mask(reference_dem)
+        ref_arr, ref_mask = get_array_and_mask(reference_dem)
 
         if inlier_mask is None:
             inlier_mask = np.ones(ref_arr.shape, dtype=bool)
@@ -1850,7 +1856,11 @@ def coregister(i: int) -> dict[str, Any] | BaseException | None:
             meta["representative_x"], meta["representative_y"] = rio.transform.xy(
                 transform_subset, representative_row, representative_col
             )
-            meta["representative_val"] = ref_subset[representative_row, representative_col]
+
+            repr_val = ref_subset[representative_row, representative_col]
+            if ~np.isfinite(repr_val):
+                repr_val = 0
+            meta["representative_val"] = repr_val
 
             # If the coreg is a pipeline, copy its metadatas to the output meta
             if hasattr(coreg, "pipeline"):
@@ -2003,12 +2013,15 @@ def subdivide_array(self, shape: tuple[int, ...]) -> NDArrayf:
         """
         if len(shape) == 3 and shape[0] == 1:  # Account for (1, row, col) shapes
             shape = (shape[1], shape[2])
-        return spatial_tools.subdivide_array(shape, count=self.subdivision)
+        return subdivide_array(shape, count=self.subdivision)
 
     def _apply_func(
         self, dem: NDArrayf, transform: rio.transform.Affine, crs: rio.crs.CRS, **kwargs: Any
     ) -> tuple[NDArrayf, rio.transform.Affine]:
 
+        if np.count_nonzero(np.isfinite(dem)) == 0:
+            return dem, transform
+
         # Other option than resample=True is not implemented for this case
         if "resample" in kwargs and kwargs["resample"] is not True:
             raise NotImplementedError()
@@ -2102,7 +2115,7 @@ def warp_dem(
     if resampling not in allowed_resampling_strs:
         raise ValueError(f"Resampling type '{resampling}' not understood. Choices: {allowed_resampling_strs}")
 
-    dem_arr, dem_mask = spatial_tools.get_array_and_mask(dem)
+    dem_arr, dem_mask = get_array_and_mask(dem)
 
     bounds, resolution = _transform_to_bounds_and_res(dem_arr.shape, transform)
 
@@ -2264,8 +2277,7 @@ def create_inlier_mask(
                 outlines = gu.Vector(shp)
             else:
                 outlines = shp
-            mask_temp = outlines.create_mask(src_dem).astype("bool")
-
+            mask_temp = outlines.create_mask(src_dem, as_array=True).reshape(np.shape(inlier_mask))
             # Append mask for given shapefile to final mask
             if inout[k] == 1:
                 inlier_mask[mask_temp] = False
@@ -2372,9 +2384,9 @@ def dem_coregistration(
 
     if isinstance(ref_dem_path, str):
         if grid == "ref":
-            ref_dem, src_dem = gu.spatial_tools.load_multiple_rasters([ref_dem_path, src_dem_path], ref_grid=0)
+            ref_dem, src_dem = gu.raster.load_multiple_rasters([ref_dem_path, src_dem_path], ref_grid=0)
         elif grid == "src":
-            ref_dem, src_dem = gu.spatial_tools.load_multiple_rasters([ref_dem_path, src_dem_path], ref_grid=1)
+            ref_dem, src_dem = gu.raster.load_multiple_rasters([ref_dem_path, src_dem_path], ref_grid=1)
     else:
         ref_dem = ref_dem_path
         src_dem = src_dem_path
diff --git a/xdem/ddem.py b/xdem/ddem.py
index f787792b..d11888a1 100644
--- a/xdem/ddem.py
+++ b/xdem/ddem.py
@@ -7,8 +7,7 @@
 import geoutils as gu
 import numpy as np
 import shapely
-from geoutils import spatial_tools
-from geoutils.georaster import RasterType
+from geoutils.raster import RasterType, get_array_and_mask
 from rasterio.crs import CRS
 from rasterio.warp import Affine
 
@@ -114,7 +113,7 @@ def from_array(
         :returns: A new dDEM instance.
         """
         return dDEM(
-            gu.georaster.Raster.from_array(data=data, transform=transform, crs=crs, nodata=nodata),
+            gu.Raster.from_array(data=data, transform=transform, crs=crs, nodata=nodata),
             start_time=start_time,
             end_time=end_time,
             error=error,
@@ -156,12 +155,12 @@ def interpolate(
             if not isinstance(mask, gu.Vector):
                 mask = gu.Vector(mask)
 
-            interpolated_ddem, nans = spatial_tools.get_array_and_mask(self.data.copy())
+            interpolated_ddem, nans = get_array_and_mask(self.data.copy())
             entries = mask.ds[mask.ds.intersects(shapely.geometry.box(*self.bounds))]
 
             ddem_mask = nans.copy().squeeze()
             for i in entries.index:
-                feature_mask = (gu.Vector(entries.loc[entries.index == i]).create_mask(self)).reshape(
+                feature_mask = (gu.Vector(entries.loc[entries.index == i]).create_mask(self, as_array=True)).reshape(
                     interpolated_ddem.shape
                 )
                 if np.count_nonzero(feature_mask) == 0:
diff --git a/xdem/dem.py b/xdem/dem.py
index 00ad79e3..dc693cd3 100644
--- a/xdem/dem.py
+++ b/xdem/dem.py
@@ -9,8 +9,8 @@
 
 import pyproj
 import rasterio as rio
-from geoutils.georaster.raster import RasterType
-from geoutils.satimg import SatelliteImage
+from geoutils import SatelliteImage
+from geoutils.raster import RasterType
 from pyproj import Transformer
 
 from xdem._typing import NDArrayf
diff --git a/xdem/demcollection.py b/xdem/demcollection.py
index d5c66026..0ec72181 100644
--- a/xdem/demcollection.py
+++ b/xdem/demcollection.py
@@ -17,10 +17,10 @@ class DEMCollection:
 
     def __init__(
         self,
-        dems: list[gu.georaster.Raster] | list[xdem.DEM],
+        dems: list[gu.Raster] | list[xdem.DEM],
         timestamps: list[datetime.datetime] | None = None,
-        outlines: gu.geovector.Vector | dict[datetime.datetime, gu.geovector.Vector] | None = None,
-        reference_dem: int | gu.georaster.Raster = 0,
+        outlines: gu.Vector | dict[datetime.datetime, gu.Vector] | None = None,
+        reference_dem: int | gu.Raster = 0,
     ):
         """
         Create a new temporal DEM collection.
@@ -55,23 +55,23 @@ def __init__(
         # The reference index changes place when sorted
         if isinstance(reference_dem, (int, np.integer)):
             self.reference_index = np.argwhere(indices == reference_dem)[0][0]
-        elif isinstance(reference_dem, gu.georaster.Raster):
+        elif isinstance(reference_dem, gu.Raster):
             self.reference_index = [i for i, dem in enumerate(self.dems) if dem is reference_dem][0]
 
         if outlines is None:
-            self.outlines: dict[np.datetime64, gu.geovector.Vector] = {}
-        elif isinstance(outlines, gu.geovector.Vector):
+            self.outlines: dict[np.datetime64, gu.Vector] = {}
+        elif isinstance(outlines, gu.Vector):
             self.outlines = {self.timestamps[self.reference_index]: outlines}
-        elif all(isinstance(value, gu.geovector.Vector) for value in outlines.values()):
+        elif all(isinstance(value, gu.Vector) for value in outlines.values()):
             self.outlines = dict(zip(np.array(list(outlines.keys())).astype("datetime64[ns]"), outlines.values()))
         else:
             raise ValueError(
                 f"Invalid format on 'outlines': {type(outlines)},"
-                " expected one of ['gu.geovector.Vector', 'dict[datetime.datetime, gu.geovector.Vector']"
+                " expected one of ['gu.Vector', 'dict[datetime.datetime, gu.Vector']"
             )
 
     @property
-    def reference_dem(self) -> gu.georaster.Raster:
+    def reference_dem(self) -> gu.Raster:
         """Get the DEM acting reference."""
         return self.dems[self.reference_index]
 
@@ -93,7 +93,7 @@ def subtract_dems(self, resampling_method: str = "cubic_spline") -> list[xdem.dD
         # Subtract every DEM that is available.
         for i, dem in enumerate(self.dems):
             # If the reference DEM is encountered, make a dDEM where dH == 0 (to keep length consistency).
-            if dem == self.reference_dem:
+            if dem.raster_equal(self.reference_dem):
                 ddem_raster = self.reference_dem.copy()
                 ddem_raster.data[:] = 0.0
                 ddem = xdem.dDEM(
@@ -154,13 +154,16 @@ def get_ddem_mask(self, ddem: xdem.dDEM, outlines_filter: str | None = None) ->
 
         # If both the start and end time outlines exist, a mask is created from their union.
         if ddem.start_time in outlines and ddem.end_time in outlines:
-            mask = np.logical_or(outlines[ddem.start_time].create_mask(ddem), outlines[ddem.end_time].create_mask(ddem))
+            mask = np.logical_or(
+                outlines[ddem.start_time].create_mask(ddem, as_array=True),
+                outlines[ddem.end_time].create_mask(ddem, as_array=True),
+            )
         # If only start time outlines exist, these should be used as a mask
         elif ddem.start_time in outlines:
-            mask = outlines[ddem.start_time].create_mask(ddem)
+            mask = outlines[ddem.start_time].create_mask(ddem, as_array=True)
         # If only one outlines file exist, use that as a mask.
         elif len(outlines) == 1:
-            mask = list(outlines.values())[0].create_mask(ddem)
+            mask = list(outlines.values())[0].create_mask(ddem, as_array=True)
         # If no fitting outlines were found, make a full true boolean mask in its stead.
         else:
             mask = np.ones(shape=ddem.data.shape, dtype=bool)
diff --git a/xdem/examples.py b/xdem/examples.py
index 1bf165a7..f0137a76 100644
--- a/xdem/examples.py
+++ b/xdem/examples.py
@@ -96,9 +96,9 @@ def process_coregistered_examples(name: str, overwrite: bool = False) -> None:
 
     # If the ddem file does not exist, create it
     if not os.path.isfile(_FILEPATHS_PROCESSED["longyearbyen_ddem"]):
-        reference_raster = gu.georaster.Raster(_FILEPATHS_DATA["longyearbyen_ref_dem"])
-        to_be_aligned_raster = gu.georaster.Raster(_FILEPATHS_DATA["longyearbyen_tba_dem"])
-        glacier_mask = gu.geovector.Vector(_FILEPATHS_DATA["longyearbyen_glacier_outlines"])
+        reference_raster = gu.Raster(_FILEPATHS_DATA["longyearbyen_ref_dem"])
+        to_be_aligned_raster = gu.Raster(_FILEPATHS_DATA["longyearbyen_tba_dem"])
+        glacier_mask = gu.Vector(_FILEPATHS_DATA["longyearbyen_glacier_outlines"])
         inlier_mask = ~glacier_mask.create_mask(reference_raster)
 
         nuth_kaab = xdem.coreg.NuthKaab()
diff --git a/xdem/fit.py b/xdem/fit.py
index f79db8c2..670cbbb0 100644
--- a/xdem/fit.py
+++ b/xdem/fit.py
@@ -10,7 +10,7 @@
 import numpy as np
 import pandas as pd
 import scipy.optimize
-from geoutils.spatial_tools import subsample_raster
+from geoutils.raster import subsample_array
 
 from xdem._typing import NDArrayf
 from xdem.spatialstats import nd_binning
@@ -279,7 +279,7 @@ def robust_polynomial_fit(
     y = y[valid_data]
 
     # Subsample data
-    subsamp = subsample_raster(x, subsample=subsample, return_indices=True, random_state=random_state)
+    subsamp = subsample_array(x, subsample=subsample, return_indices=True, random_state=random_state)
     x = x[subsamp]
     y = y[subsamp]
 
@@ -454,7 +454,7 @@ def wrapper_cost_sumofsin(p: NDArrayf, x: NDArrayf, y: NDArrayf) -> float:
         init_x = np.array([np.round(ini, 5) for ini in init_results.x])
 
         # Subsample the final raster
-        subsamp = subsample_raster(x, subsample=subsample, return_indices=True, random_state=random_state)
+        subsamp = subsample_array(x, subsample=subsample, return_indices=True, random_state=random_state)
         x = x[subsamp]
         y = y[subsamp]
 
diff --git a/xdem/misc.py b/xdem/misc.py
index 12c89569..43cf2aac 100644
--- a/xdem/misc.py
+++ b/xdem/misc.py
@@ -1,6 +1,7 @@
 """Small functions for testing, examples, and other miscellaneous uses."""
 from __future__ import annotations
 
+import copy
 import functools
 import warnings
 from typing import Any, Callable
@@ -125,21 +126,30 @@ def new_func(*args: Any, **kwargs: Any) -> Any:
     return deprecator_func
 
 
-def diff_environment_yml(fn_env: str, fn_devenv: str, print_dep: str = "both") -> None:
+def diff_environment_yml(
+    fn_env: str | dict[str, Any], fn_devenv: str | dict[str, Any], print_dep: str = "both", input_dict: bool = False
+) -> None:
     """
     Compute the difference between environment.yml and dev-environment.yml for setup of continuous integration,
     while checking that all the dependencies listed in environment.yml are also in dev-environment.yml
     :param fn_env: Filename path to environment.yml
     :param fn_devenv: Filename path to dev-environment.yml
     :param print_dep: Whether to print conda differences "conda", pip differences "pip" or both.
+    :param input_dict: Whether to consider the input as a dict (for testing purposes).
     """
 
     if not _has_yaml:
         raise ValueError("Test dependency needed. Install 'pyyaml'")
 
-    # Load the yml as dictionaries
-    yaml_env = yaml.safe_load(open(fn_env))
-    yaml_devenv = yaml.safe_load(open(fn_devenv))
+    if not input_dict:
+        # Load the yml as dictionaries
+        yaml_env = yaml.safe_load(open(fn_env))  # type: ignore
+        yaml_devenv = yaml.safe_load(open(fn_devenv))  # type: ignore
+    else:
+        # We need a copy as we'll pop things out and don't want to affect input
+        # dict.copy() is shallow and does not work with embedded list in dicts (as is the case here)
+        yaml_env = copy.deepcopy(fn_env)
+        yaml_devenv = copy.deepcopy(fn_devenv)
 
     # Extract the dependencies values
     conda_dep_env = yaml_env["dependencies"]
@@ -149,6 +159,10 @@ def diff_environment_yml(fn_env: str, fn_devenv: str, print_dep: str = "both") -
     if isinstance(conda_dep_devenv[-1], dict):
         pip_dep_devenv = conda_dep_devenv.pop()["pip"]
 
+        # Remove the package's self install for devs via pip, if it exists
+        if "-e ./" in pip_dep_devenv:
+            pip_dep_devenv.remove("-e ./")
+
         # Check if there is a pip dependency in the normal env as well, if yes pop it also
         if isinstance(conda_dep_env[-1], dict):
             pip_dep_env = conda_dep_env.pop()["pip"]
@@ -166,7 +180,7 @@ def diff_environment_yml(fn_env: str, fn_devenv: str, print_dep: str = "both") -
 
         # If there is no pip dependency in env, all the ones of dev-env need to be added during CI
         else:
-            diff_pip_dep = list(pip_dep_devenv["pip"])
+            diff_pip_dep = pip_dep_devenv
 
     # If there is no pip dependency, we ignore this step
     else:
diff --git a/xdem/spatialstats.py b/xdem/spatialstats.py
index 1d5b28bf..533e84f0 100644
--- a/xdem/spatialstats.py
+++ b/xdem/spatialstats.py
@@ -15,9 +15,14 @@
 import numba
 import numpy as np
 import pandas as pd
-from geoutils.georaster import Raster, RasterType
-from geoutils.geovector import Vector, VectorType
-from geoutils.spatial_tools import get_array_and_mask, subsample_raster
+from geoutils.raster import (
+    Mask,
+    Raster,
+    RasterType,
+    get_array_and_mask,
+    subsample_array,
+)
+from geoutils.vector import Vector, VectorType
 from numba import jit
 from numpy.typing import ArrayLike
 from scipy import integrate
@@ -35,20 +40,22 @@
     import skgstat as skg
 
 
-def nmad(data: NDArrayf, nfact: float = 1.4826) -> np.floating[Any]:
+def nmad(data: NDArrayf | RasterType, nfact: float = 1.4826) -> np.floating[Any]:
     """
     Calculate the normalized median absolute deviation (NMAD) of an array.
     Default scaling factor is 1.4826 to scale the median absolute deviation (MAD) to the dispersion of a normal
     distribution (see https://en.wikipedia.org/wiki/Median_absolute_deviation#Relation_to_standard_deviation, and
     e.g. Höhle and Höhle (2009), http://dx.doi.org/10.1016/j.isprsjprs.2009.02.003)
 
-    :param data: Input data
+    :param data: Input array or raster
     :param nfact: Normalization factor for the data
 
     :returns nmad: (normalized) median absolute deviation of data.
     """
     if isinstance(data, np.ma.masked_array):
         data_arr = get_array_and_mask(data, check_shape=False)[0]
+    elif isinstance(data, Raster):
+        data_arr = data
     else:
         data_arr = np.asarray(data)
     return nfact * np.nanmedian(np.abs(data_arr - np.nanmedian(data_arr)))
@@ -444,8 +451,8 @@ def estimate_model_heteroscedasticity(
 @overload
 def _preprocess_values_with_mask_to_array(  # type: ignore
     values: list[NDArrayf | RasterType],
-    include_mask: NDArrayf | VectorType | gpd.GeoDataFrame = None,
-    exclude_mask: NDArrayf | VectorType | gpd.GeoDataFrame = None,
+    include_mask: NDArrayf | Mask | VectorType | gpd.GeoDataFrame = None,
+    exclude_mask: NDArrayf | Mask | VectorType | gpd.GeoDataFrame = None,
     gsd: float | None = None,
     preserve_shape: bool = True,
 ) -> tuple[list[NDArrayf], float]:
@@ -455,8 +462,8 @@ def _preprocess_values_with_mask_to_array(  # type: ignore
 @overload
 def _preprocess_values_with_mask_to_array(
     values: NDArrayf | RasterType,
-    include_mask: NDArrayf | VectorType | gpd.GeoDataFrame = None,
-    exclude_mask: NDArrayf | VectorType | gpd.GeoDataFrame = None,
+    include_mask: NDArrayf | Mask | VectorType | gpd.GeoDataFrame = None,
+    exclude_mask: NDArrayf | Mask | VectorType | gpd.GeoDataFrame = None,
     gsd: float | None = None,
     preserve_shape: bool = True,
 ) -> tuple[NDArrayf, float]:
@@ -465,8 +472,8 @@ def _preprocess_values_with_mask_to_array(
 
 def _preprocess_values_with_mask_to_array(
     values: list[NDArrayf | RasterType] | NDArrayf | RasterType,
-    include_mask: NDArrayf | VectorType | gpd.GeoDataFrame = None,
-    exclude_mask: NDArrayf | VectorType | gpd.GeoDataFrame = None,
+    include_mask: NDArrayf | Mask | VectorType | gpd.GeoDataFrame = None,
+    exclude_mask: NDArrayf | Mask | VectorType | gpd.GeoDataFrame = None,
     gsd: float | None = None,
     preserve_shape: bool = True,
 ) -> tuple[list[NDArrayf] | NDArrayf, float]:
@@ -493,10 +500,10 @@ def _preprocess_values_with_mask_to_array(
     ):
         raise ValueError("The values must be a Raster or NumPy array, or a list of those.")
     # Masks need to be an array, Vector or GeoPandas dataframe
-    if include_mask is not None and not isinstance(include_mask, (np.ndarray, Vector, gpd.GeoDataFrame)):
-        raise ValueError("The stable mask must be a Vector, GeoDataFrame or NumPy array.")
-    if exclude_mask is not None and not isinstance(exclude_mask, (np.ndarray, Vector, gpd.GeoDataFrame)):
-        raise ValueError("The unstable mask must be a Vector, GeoDataFrame or NumPy array.")
+    if include_mask is not None and not isinstance(include_mask, (np.ndarray, Vector, Mask, gpd.GeoDataFrame)):
+        raise ValueError("The stable mask must be a Vector, Mask, GeoDataFrame or NumPy array.")
+    if exclude_mask is not None and not isinstance(exclude_mask, (np.ndarray, Vector, Mask, gpd.GeoDataFrame)):
+        raise ValueError("The unstable mask must be a Vector, Mask, GeoDataFrame or NumPy array.")
 
     # Check that input stable mask can only be a georeferenced vector if the proxy values are a Raster to project onto
     if isinstance(values, list):
@@ -543,7 +550,10 @@ def _preprocess_values_with_mask_to_array(
             stable_vector = include_mask
 
         # Create the mask
-        include_mask_arr = stable_vector.create_mask(first_raster)
+        include_mask_arr = stable_vector.create_mask(first_raster, as_array=True)
+    # If the mask is a Mask
+    elif isinstance(include_mask, Mask):
+        include_mask_arr = include_mask.data.filled(False)
     # If the mask is already an array, just pass it
     else:
         include_mask_arr = include_mask
@@ -560,8 +570,11 @@ def _preprocess_values_with_mask_to_array(
             unstable_vector = exclude_mask
 
         # Create the mask
-        exclude_mask_arr = unstable_vector.create_mask(first_raster)
+        exclude_mask_arr = unstable_vector.create_mask(first_raster, as_array=True)
     # If the mask is already an array, just pass it
+    # If the mask is a Mask
+    elif isinstance(exclude_mask, Mask):
+        exclude_mask_arr = exclude_mask.data.filled(False)
     else:
         exclude_mask_arr = exclude_mask
 
@@ -588,8 +601,8 @@ def _preprocess_values_with_mask_to_array(
 def infer_heteroscedasticity_from_stable(
     dvalues: NDArrayf,
     list_var: list[NDArrayf | RasterType],
-    stable_mask: NDArrayf | VectorType | gpd.GeoDataFrame = None,
-    unstable_mask: NDArrayf | VectorType | gpd.GeoDataFrame = None,
+    stable_mask: NDArrayf | Mask | VectorType | gpd.GeoDataFrame = None,
+    unstable_mask: NDArrayf | Mask | VectorType | gpd.GeoDataFrame = None,
     list_var_names: list[str] = None,
     spread_statistic: Callable[[NDArrayf], np.floating[Any]] = nmad,
     list_var_bins: int | tuple[int, ...] | tuple[NDArrayf] | None = None,
@@ -603,8 +616,8 @@ def infer_heteroscedasticity_from_stable(
 def infer_heteroscedasticity_from_stable(
     dvalues: RasterType,
     list_var: list[NDArrayf | RasterType],
-    stable_mask: NDArrayf | VectorType | gpd.GeoDataFrame = None,
-    unstable_mask: NDArrayf | VectorType | gpd.GeoDataFrame = None,
+    stable_mask: NDArrayf | Mask | VectorType | gpd.GeoDataFrame = None,
+    unstable_mask: NDArrayf | Mask | VectorType | gpd.GeoDataFrame = None,
     list_var_names: list[str] = None,
     spread_statistic: Callable[[NDArrayf], np.floating[Any]] = nmad,
     list_var_bins: int | tuple[int, ...] | tuple[NDArrayf] | None = None,
@@ -617,8 +630,8 @@ def infer_heteroscedasticity_from_stable(
 def infer_heteroscedasticity_from_stable(
     dvalues: NDArrayf | RasterType,
     list_var: list[NDArrayf | RasterType],
-    stable_mask: NDArrayf | VectorType | gpd.GeoDataFrame = None,
-    unstable_mask: NDArrayf | VectorType | gpd.GeoDataFrame = None,
+    stable_mask: NDArrayf | Mask | VectorType | gpd.GeoDataFrame = None,
+    unstable_mask: NDArrayf | Mask | VectorType | gpd.GeoDataFrame = None,
     list_var_names: list[str] = None,
     spread_statistic: Callable[[NDArrayf], np.floating[Any]] = nmad,
     list_var_bins: int | tuple[int, ...] | tuple[NDArrayf] | None = None,
@@ -810,7 +823,7 @@ def _subsample_wrapper(
         values_sp = values
         coords_sp = coords
 
-    index = subsample_raster(values_sp, subsample=subsample, return_indices=True, random_state=rnd)
+    index = subsample_array(values_sp, subsample=subsample, return_indices=True, random_state=rnd)
     values_sub = values_sp[index[0]]
     coords_sub = coords_sp[index[0], :]
 
@@ -1696,8 +1709,8 @@ def estimate_model_spatial_correlation(
 def infer_spatial_correlation_from_stable(
     dvalues: NDArrayf | RasterType,
     list_models: list[str | Callable[[NDArrayf, float, float], NDArrayf]],
-    stable_mask: NDArrayf | VectorType | gpd.GeoDataFrame = None,
-    unstable_mask: NDArrayf | VectorType | gpd.GeoDataFrame = None,
+    stable_mask: NDArrayf | Mask | VectorType | gpd.GeoDataFrame = None,
+    unstable_mask: NDArrayf | Mask | VectorType | gpd.GeoDataFrame = None,
     errors: NDArrayf | RasterType = None,
     estimator: str = "dowd",
     gsd: float = None,
@@ -2212,7 +2225,7 @@ def number_effective_samples(
         if isinstance(rasterize_resolution, (float, int, np.floating, np.integer)):
 
             # We only need relative mask and coordinates, not absolute
-            mask = V.create_mask(xres=rasterize_resolution)
+            mask = V.create_mask(xres=rasterize_resolution, as_array=True)
             x = rasterize_resolution * np.arange(0, mask.shape[0])
             y = rasterize_resolution * np.arange(0, mask.shape[1])
             coords = np.array(np.meshgrid(y, x))
@@ -2221,7 +2234,7 @@ def number_effective_samples(
         elif isinstance(rasterize_resolution, Raster):
 
             # With a Raster we can get the coordinates directly
-            mask = V.create_mask(rst=rasterize_resolution).squeeze()
+            mask = V.create_mask(rst=rasterize_resolution, as_array=True).squeeze()
             coords = np.array(rasterize_resolution.coords())
             coords_on_mask = coords[:, mask].T
 
@@ -2286,7 +2299,7 @@ def spatial_error_propagation(
                 area_vector = Vector(area)
             else:
                 area_vector = area
-            area_mask = area_vector.create_mask(errors).squeeze()
+            area_mask = area_vector.create_mask(errors, as_array=True).squeeze()
 
             average_spread = np.nanmean(errors_arr[area_mask])
 
diff --git a/xdem/terrain.py b/xdem/terrain.py
index e9def76b..f2884aa0 100644
--- a/xdem/terrain.py
+++ b/xdem/terrain.py
@@ -7,7 +7,7 @@
 import geoutils as gu
 import numba
 import numpy as np
-from geoutils.georaster import Raster, RasterType
+from geoutils.raster import Raster, RasterType
 
 from xdem._typing import MArrayf, NDArrayf
 
@@ -300,7 +300,7 @@ def get_quadric_coefficients(
     :returns: An array of coefficients for each pixel of shape (9, row, col).
     """
     # This function only formats and validates the inputs. For the true functionality, see _get_quadric_coefficients()
-    dem_arr = gu.spatial_tools.get_array_and_mask(dem)[0]
+    dem_arr = gu.raster.get_array_and_mask(dem)[0]
 
     if len(dem_arr.shape) != 2:
         raise ValueError(
@@ -572,7 +572,7 @@ def get_windowed_indexes(
     :returns: An array of coefficients for each pixel of shape (5, row, col).
     """
     # This function only formats and validates the inputs. For the true functionality, see _get_quadric_coefficients()
-    dem_arr = gu.spatial_tools.get_array_and_mask(dem)[0]
+    dem_arr = gu.raster.get_array_and_mask(dem)[0]
 
     if len(dem_arr.shape) != 2:
         raise ValueError(
@@ -873,7 +873,7 @@ def get_terrain_attribute(
     make_fractal_roughness = "fractal_roughness" in attribute
 
     # Get array of DEM
-    dem_arr = gu.spatial_tools.get_array_and_mask(dem)[0]
+    dem_arr = gu.raster.get_array_and_mask(dem)[0]
 
     if make_surface_fit:
         if not isinstance(resolution, Sized):
diff --git a/xdem/volume.py b/xdem/volume.py
index 267ae6fb..21784309 100644
--- a/xdem/volume.py
+++ b/xdem/volume.py
@@ -10,7 +10,7 @@
 import pandas as pd
 import rasterio.fill
 import scipy.interpolate
-from geoutils import spatial_tools
+from geoutils.raster import RasterType, get_array_and_mask, get_mask, get_valid_extent
 from tqdm import tqdm
 
 import xdem
@@ -42,10 +42,10 @@ def hypsometric_binning(
     assert ddem.shape == ref_dem.shape
 
     # Convert ddem mask into NaN
-    ddem, _ = spatial_tools.get_array_and_mask(ddem)
+    ddem, _ = get_array_and_mask(ddem)
 
     # Extract only the valid values, i.e. valid in ref_dem
-    valid_mask = ~spatial_tools.get_mask(ref_dem)
+    valid_mask = ~get_mask(ref_dem)
     ddem = np.array(ddem[valid_mask])
     ref_dem = np.array(ref_dem.squeeze()[valid_mask])
 
@@ -290,7 +290,7 @@ def linear_interpolation(
     :returns: A filled array with no NaNs
     """
     # Create a mask for where nans exist
-    nan_mask = spatial_tools.get_mask(array)
+    nan_mask = get_mask(array)
 
     interpolated_array = rasterio.fill.fillnodata(
         array.copy(), mask=(~nan_mask).astype("uint8"), max_search_distance=max_search_distance
@@ -338,10 +338,10 @@ def hypsometric_interpolation(
     :param mask: A mask to delineate the area that will be interpolated (True means hypsometric will be used).
     """
     # Get ddem array with invalid pixels converted to NaN and mask of invalid pixels
-    ddem, ddem_mask = spatial_tools.get_array_and_mask(voided_ddem)
+    ddem, ddem_mask = get_array_and_mask(voided_ddem)
 
     # Get ref_dem array with invalid pixels converted to NaN and mask of invalid pixels
-    dem, dem_mask = spatial_tools.get_array_and_mask(ref_dem)
+    dem, dem_mask = get_array_and_mask(ref_dem)
 
     # Make sure the mask does not have e.g. the shape (1, height, width)
     mask = mask.squeeze()
@@ -413,10 +413,10 @@ def local_hypsometric_interpolation(
     assert voided_ddem.shape == ref_dem.shape == mask.shape
 
     # Get ddem array with invalid pixels converted to NaN and mask of invalid pixels
-    ddem, ddem_mask = spatial_tools.get_array_and_mask(voided_ddem)
+    ddem, ddem_mask = get_array_and_mask(voided_ddem)
 
     # Get ref_dem array with invalid pixels converted to NaN and mask of invalid pixels
-    dem, dem_mask = spatial_tools.get_array_and_mask(ref_dem)
+    dem, dem_mask = get_array_and_mask(ref_dem)
 
     # A mask of inlier values: The union of the mask and the inverted exclusion masks of both rasters.
     inlier_mask = (mask != 0) & (~ddem_mask & ~dem_mask)
@@ -482,7 +482,7 @@ def local_hypsometric_interpolation(
         if plot:
             local_ddem = np.where(local_inlier_mask, ddem, np.nan)
             vmax = max(np.abs(np.nanpercentile(local_ddem, [2, 98])))
-            rowmin, rowmax, colmin, colmax = spatial_tools.get_valid_extent(mask == index)
+            rowmin, rowmax, colmin, colmax = get_valid_extent(mask == index)
 
             plt.figure(figsize=(12, 8))
             plt.subplot(121)
@@ -534,9 +534,9 @@ def local_hypsometric_interpolation(
 
 
 def get_regional_hypsometric_signal(
-    ddem: NDArrayf | MArrayf,
-    ref_dem: NDArrayf | MArrayf,
-    glacier_index_map: NDArrayf,
+    ddem: NDArrayf | MArrayf | RasterType,
+    ref_dem: NDArrayf | MArrayf | RasterType,
+    glacier_index_map: NDArrayf | RasterType,
     n_bins: int = 20,
     verbose: bool = False,
     min_coverage: float = 0.05,
@@ -554,8 +554,9 @@ def get_regional_hypsometric_signal(
     :returns: A DataFrame of bin statistics, scaled by elevation and elevation change.
     """
     # Extract the array and mask representations of the arrays.
-    ddem_arr, ddem_mask = spatial_tools.get_array_and_mask(ddem.squeeze())
-    ref_arr, ref_mask = spatial_tools.get_array_and_mask(ref_dem.squeeze())
+    ddem_arr, ddem_mask = get_array_and_mask(ddem)
+    ref_arr, ref_mask = get_array_and_mask(ref_dem)
+    glacier_index_map, _ = get_array_and_mask(glacier_index_map)
 
     # The reference DEM should be void free
     assert np.count_nonzero(ref_mask) == 0, "Reference DEM has voids"
@@ -631,9 +632,9 @@ def get_regional_hypsometric_signal(
 
 
 def norm_regional_hypsometric_interpolation(
-    voided_ddem: NDArrayf | MArrayf,
-    ref_dem: NDArrayf | MArrayf,
-    glacier_index_map: NDArrayf,
+    voided_ddem: NDArrayf | MArrayf | RasterType,
+    ref_dem: NDArrayf | MArrayf | RasterType,
+    glacier_index_map: NDArrayf | RasterType,
     min_coverage: float = 0.1,
     regional_signal: pd.DataFrame | None = None,
     verbose: bool = False,
@@ -660,8 +661,9 @@ def norm_regional_hypsometric_interpolation(
     :returns: A dDEM where glacier's that fit the min_coverage criterion are interpolated.
     """
     # Extract the array and nan parts of the inputs.
-    ddem_arr, ddem_nans = spatial_tools.get_array_and_mask(voided_ddem)
-    ref_arr, ref_nans = spatial_tools.get_array_and_mask(ref_dem)
+    ddem_arr, ddem_nans = get_array_and_mask(voided_ddem)
+    ref_arr, ref_nans = get_array_and_mask(ref_dem)
+    glacier_index_map, _ = get_array_and_mask(glacier_index_map)
 
     # The reference DEM should be void free
     assert np.count_nonzero(ref_nans) == 0, "Reference DEM has voids"