From 2ae17335e91e307b5d74f83f1a85a84b4724fedb Mon Sep 17 00:00:00 2001
From: Romain Hugonnet <romain.hugonnet@gmail.com>
Date: Wed, 13 Nov 2024 21:05:50 -0900
Subject: [PATCH 01/10] Update with GeoUtils restructuration

---
 tests/test_coreg/test_affine.py |  4 ++--
 xdem/coreg/affine.py            |  2 +-
 xdem/coreg/base.py              | 25 ++++++++++++-------------
 xdem/ddem.py                    |  3 ++-
 xdem/spatialstats.py            | 11 +++--------
 xdem/volume.py                  |  4 ++--
 6 files changed, 22 insertions(+), 27 deletions(-)

diff --git a/tests/test_coreg/test_affine.py b/tests/test_coreg/test_affine.py
index 9e82d556..86f1117d 100644
--- a/tests/test_coreg/test_affine.py
+++ b/tests/test_coreg/test_affine.py
@@ -12,7 +12,7 @@
 from geoutils import Raster, Vector
 from geoutils._typing import NDArrayNum
 from geoutils.raster import RasterType
-from geoutils.raster.raster import _shift_transform
+from geoutils.raster.geotransformations import _translate
 from scipy.ndimage import binary_dilation
 
 from xdem import coreg, examples
@@ -128,7 +128,7 @@ def test_reproject_horizontal_shift_samecrs__gdal(self, xoff_yoff: tuple[float,
 
         # Reproject with SciPy
         xoff, yoff = xoff_yoff
-        dst_transform = _shift_transform(transform=ref.transform, xoff=xoff, yoff=yoff, distance_unit="georeferenced")
+        dst_transform = _translate(transform=ref.transform, xoff=xoff, yoff=yoff, distance_unit="georeferenced")
         output = _reproject_horizontal_shift_samecrs(
             raster_arr=ref.data, src_transform=ref.transform, dst_transform=dst_transform
         )
diff --git a/xdem/coreg/affine.py b/xdem/coreg/affine.py
index 69090979..31bca945 100644
--- a/xdem/coreg/affine.py
+++ b/xdem/coreg/affine.py
@@ -18,8 +18,8 @@
 import numpy as np
 import rasterio as rio
 import scipy.optimize
+from geoutils.interface.interpolate import _interp_points
 from geoutils.raster.georeferencing import _bounds, _coords, _res
-from geoutils.raster.interpolate import _interp_points
 from tqdm import trange
 
 from xdem._typing import NDArrayb, NDArrayf
diff --git a/xdem/coreg/base.py b/xdem/coreg/base.py
index 7ebceb1a..9d26f75d 100644
--- a/xdem/coreg/base.py
+++ b/xdem/coreg/base.py
@@ -32,17 +32,17 @@
 import scipy.optimize
 import skimage.transform
 from geoutils._typing import Number
-from geoutils.misc import resampling_method_from_str
-from geoutils.raster import (
-    Mask,
-    RasterType,
-    get_array_and_mask,
-    raster,
-    subdivide_array,
+from geoutils.interface.gridding import _grid_pointcloud
+from geoutils.interface.interpolate import _interp_points
+from geoutils.raster import Mask, RasterType, raster, subdivide_array
+from geoutils.raster.array import get_array_and_mask
+from geoutils.raster.georeferencing import (
+    _bounds,
+    _cast_pixel_interpretation,
+    _coords,
+    _res,
 )
-from geoutils.raster.georeferencing import _bounds, _coords, _res
-from geoutils.raster.interpolate import _interp_points
-from geoutils.raster.raster import _cast_pixel_interpretation, _shift_transform
+from geoutils.raster.geotransformations import _resampling_method_from_str, _translate
 from tqdm import tqdm
 
 from xdem._typing import MArrayf, NDArrayb, NDArrayf
@@ -523,7 +523,7 @@ def _postprocess_coreg_apply(
     """
 
     # Define resampling
-    resampling = resampling if isinstance(resampling, rio.warp.Resampling) else resampling_method_from_str(resampling)
+    resampling = resampling if isinstance(resampling, rio.warp.Resampling) else _resampling_method_from_str(resampling)
 
     # Distribute between raster and point apply methods
     if isinstance(applied_elev, np.ndarray):
@@ -1244,7 +1244,7 @@ def _apply_matrix_rst(
 
     # 2/ Check if the matrix contains only translations, in that case only shift the DEM only by translation
     if np.array_equal(shift_only_matrix, matrix) and force_regrid_method is None:
-        new_transform = _shift_transform(transform, xoff=matrix[0, 3], yoff=matrix[1, 3])
+        new_transform = _translate(transform, xoff=matrix[0, 3], yoff=matrix[1, 3])
         return dem + matrix[2, 3], new_transform
 
     # 3/ If matrix contains only small rotations (less than 20 degrees), use the fast iterative reprojection
@@ -1260,7 +1260,6 @@ def _apply_matrix_rst(
     dem_rst = gu.Raster.from_array(dem, transform=transform, crs=None, nodata=99999)
     epc = dem_rst.to_pointcloud(data_column_name="z").ds
     trans_epc = _apply_matrix_pts(epc, matrix=matrix, centroid=centroid)
-    from geoutils.pointcloud import _grid_pointcloud
 
     new_dem = _grid_pointcloud(
         trans_epc, grid_coords=dem_rst.coords(grid=False), data_column_name="z", resampling=resampling
diff --git a/xdem/ddem.py b/xdem/ddem.py
index 487fe903..0a2ce95e 100644
--- a/xdem/ddem.py
+++ b/xdem/ddem.py
@@ -9,7 +9,8 @@
 import pyogrio
 import rasterio as rio
 import shapely
-from geoutils.raster import Raster, RasterType, get_array_and_mask
+from geoutils.raster import Raster, RasterType
+from geoutils.raster.array import get_array_and_mask
 from rasterio.crs import CRS
 from rasterio.warp import Affine
 
diff --git a/xdem/spatialstats.py b/xdem/spatialstats.py
index 1898fb62..0602cd1e 100644
--- a/xdem/spatialstats.py
+++ b/xdem/spatialstats.py
@@ -16,14 +16,9 @@
 import numba
 import numpy as np
 import pandas as pd
-from geoutils.raster import (
-    Mask,
-    Raster,
-    RasterType,
-    get_array_and_mask,
-    subsample_array,
-)
-from geoutils.vector import Vector, VectorType
+from geoutils.raster import Mask, Raster, RasterType, subsample_array
+from geoutils.raster.array import get_array_and_mask
+from geoutils.vector.vector import Vector, VectorType
 from numpy.typing import ArrayLike
 from scipy import integrate
 from scipy.interpolate import RegularGridInterpolator, griddata
diff --git a/xdem/volume.py b/xdem/volume.py
index 65d09447..bf9ae67b 100644
--- a/xdem/volume.py
+++ b/xdem/volume.py
@@ -10,8 +10,8 @@
 import pandas as pd
 import rasterio.fill
 import scipy.interpolate
-from geoutils.raster import (
-    RasterType,
+from geoutils.raster import RasterType
+from geoutils.raster.array import (
     get_array_and_mask,
     get_mask_from_array,
     get_valid_extent,

From d5699e7b8f7e616ff9660e0c69df94093f56e933 Mon Sep 17 00:00:00 2001
From: Romain Hugonnet <romain.hugonnet@gmail.com>
Date: Wed, 13 Nov 2024 21:08:15 -0900
Subject: [PATCH 02/10] Bump GeoUtils version number in envs

---
 dev-environment.yml | 2 +-
 environment.yml     | 2 +-
 requirements.txt    | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/dev-environment.yml b/dev-environment.yml
index fb682146..78004d13 100644
--- a/dev-environment.yml
+++ b/dev-environment.yml
@@ -13,7 +13,7 @@ dependencies:
   - tqdm
   - scikit-image=0.*
   - scikit-gstat>=1.0.18,<1.1
-  - geoutils=0.1.10
+  - geoutils=0.1.11
 
   # Development-specific, to mirror manually in setup.cfg [options.extras_require].
   - pip
diff --git a/environment.yml b/environment.yml
index 70613281..55e391e7 100644
--- a/environment.yml
+++ b/environment.yml
@@ -13,7 +13,7 @@ dependencies:
   - tqdm
   - scikit-image=0.*
   - scikit-gstat>=1.0.18,<1.1
-  - geoutils=0.1.10
+  - geoutils=0.1.11
   - pip
 
    # To run CI against latest GeoUtils
diff --git a/requirements.txt b/requirements.txt
index 0d5001f8..d8f6ecfd 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -11,5 +11,5 @@ scipy==1.*
 tqdm
 scikit-image==0.*
 scikit-gstat>=1.0.18,<1.1
-geoutils==0.1.10
+geoutils==0.1.11
 pip

From 9f0c72036be51fbf5555be656798544fe94a7f01 Mon Sep 17 00:00:00 2001
From: Romain Hugonnet <romain.hugonnet@gmail.com>
Date: Sun, 17 Nov 2024 16:13:51 -0900
Subject: [PATCH 03/10] Update with removal SatelliteImage

---
 doc/source/vertical_ref.md |  6 +++---
 xdem/dem.py                | 22 +++++++++++++++-------
 2 files changed, 18 insertions(+), 10 deletions(-)

diff --git a/doc/source/vertical_ref.md b/doc/source/vertical_ref.md
index 3c49f2de..460e3a61 100644
--- a/doc/source/vertical_ref.md
+++ b/doc/source/vertical_ref.md
@@ -130,11 +130,11 @@ import os
 os.remove("mydem_with3dcrs.tif")
 ```
 
-2. **From the {attr}`~xdem.DEM.product` name of the DEM**, if parsed from the filename by {class}`geoutils.SatelliteImage`.
+2. **From the {attr}`~xdem.DEM.product` name of the DEM**, if parsed from the filename by the ``parse_sensor_metadata`` of {class}`geoutils.Raster`.
 
 
 ```{seealso}
-The {class}`~geoutils.SatelliteImage` parent class that parses the product metadata is described in [a dedicated section of GeoUtils' documentation](https://geoutils.readthedocs.io/en/latest/satimg_class.html).
+The {class}`~geoutils.Raster` parent class can parse sensor metadata, see its [documentation page](https://geoutils.readthedocs.io/en/stable/core_satimg.html).
 ```
 
 ```{code-cell} ipython3
@@ -151,7 +151,7 @@ dem.save("SETSM_WV03_20151101_104001001327F500_104001001312DE00_seg2_2m_v3.0_dem
 ```
 
 ```{code-cell} ipython3
-# Open an ArcticDEM strip file, recognized as an ArcticDEM product by SatelliteImage
+# Open an ArcticDEM strip file, recognized as an ArcticDEM product
 dem = xdem.DEM("SETSM_WV03_20151101_104001001327F500_104001001312DE00_seg2_2m_v3.0_dem.tif")
 # The vertical CRS is set as "Ellipsoid" from knowing that is it an ArcticDEM product
 dem.vcrs
diff --git a/xdem/dem.py b/xdem/dem.py
index 7d6533d8..561a782a 100644
--- a/xdem/dem.py
+++ b/xdem/dem.py
@@ -1,4 +1,4 @@
-"""DEM class and functions."""
+"""This module defines the DEM class."""
 from __future__ import annotations
 
 import pathlib
@@ -9,7 +9,7 @@
 import numpy as np
 import rasterio as rio
 from affine import Affine
-from geoutils import SatelliteImage
+from geoutils import Raster
 from geoutils.raster import Mask, RasterType
 from pyproj import CRS
 from pyproj.crs import CompoundCRS, VerticalCRS
@@ -35,7 +35,7 @@
 dem_attrs = ["_vcrs", "_vcrs_name", "_vcrs_grid"]
 
 
-class DEM(SatelliteImage):  # type: ignore
+class DEM(Raster):  # type: ignore
     """
     The digital elevation model.
 
@@ -76,8 +76,11 @@ def __init__(
         | pathlib.Path
         | int
         | None = None,
+        load_data: bool = False,
+        parse_sensor_metadata: bool = False,
         silent: bool = True,
-        **kwargs: Any,
+        downsample: int = 1,
+        nodata: int | float | None = None,
     ) -> None:
         """
         Instantiate a digital elevation model.
@@ -90,7 +93,11 @@ def __init__(
         :param filename_or_dataset: The filename of the dataset.
         :param vcrs: Vertical coordinate reference system either as a name ("WGS84", "EGM08", "EGM96"),
             an EPSG code or pyproj.crs.VerticalCRS, or a path to a PROJ grid file (https://github.com/OSGeo/PROJ-data).
+        :param load_data: Whether to load the array during instantiation. Default is False.
+        :param parse_sensor_metadata: Whether to parse sensor metadata from filename and similarly-named metadata files.
         :param silent: Whether to display vertical reference parsing.
+        :param downsample: Downsample the array once loaded by a round factor. Default is no downsampling.
+        :param nodata: Nodata value to be used (overwrites the metadata). Default reads from metadata.
         """
 
         self.data: NDArrayf
@@ -107,7 +114,8 @@ def __init__(
         else:
             with warnings.catch_warnings():
                 warnings.filterwarnings("ignore", message="Parse metadata from file not implemented")
-                super().__init__(filename_or_dataset, silent=silent, **kwargs)
+                super().__init__(filename_or_dataset, load_data=load_data, parse_sensor_metadata=parse_sensor_metadata,
+                                 silent=silent, downsample=downsample, nodata=nodata)
 
         # Ensure DEM has only one band: self.bands can be None when data is not loaded through the Raster class
         if self.bands is not None and len(self.bands) > 1:
@@ -134,8 +142,8 @@ def __init__(
                 vcrs = vcrs_from_crs
 
         # If no vertical CRS was provided by the user or defined in the CRS
-        if vcrs is None:
-            vcrs = _parse_vcrs_name_from_product(self.product)
+        if vcrs is None and "product" in self.tags:
+            vcrs = _parse_vcrs_name_from_product(self.tags["product"])
 
         # If a vertical reference was parsed or provided by user
         if vcrs is not None:

From c177d8bde6ee816c7fc55fb6e8db0dba5b48fe65 Mon Sep 17 00:00:00 2001
From: Romain Hugonnet <romain.hugonnet@gmail.com>
Date: Sun, 17 Nov 2024 16:14:31 -0900
Subject: [PATCH 04/10] Linting

---
 xdem/dem.py | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/xdem/dem.py b/xdem/dem.py
index 561a782a..9cbe0987 100644
--- a/xdem/dem.py
+++ b/xdem/dem.py
@@ -114,8 +114,14 @@ def __init__(
         else:
             with warnings.catch_warnings():
                 warnings.filterwarnings("ignore", message="Parse metadata from file not implemented")
-                super().__init__(filename_or_dataset, load_data=load_data, parse_sensor_metadata=parse_sensor_metadata,
-                                 silent=silent, downsample=downsample, nodata=nodata)
+                super().__init__(
+                    filename_or_dataset,
+                    load_data=load_data,
+                    parse_sensor_metadata=parse_sensor_metadata,
+                    silent=silent,
+                    downsample=downsample,
+                    nodata=nodata,
+                )
 
         # Ensure DEM has only one band: self.bands can be None when data is not loaded through the Raster class
         if self.bands is not None and len(self.bands) > 1:
@@ -200,7 +206,7 @@ def from_array(
         :returns: DEM created from the provided array and georeferencing.
         """
         # We first apply the from_array of the parent class
-        rast = SatelliteImage.from_array(
+        rast = Raster.from_array(
             data=data,
             transform=transform,
             crs=crs,

From cf3da605dd19942a9a4cb56246ca7919a203d934 Mon Sep 17 00:00:00 2001
From: Romain Hugonnet <romain.hugonnet@gmail.com>
Date: Sun, 17 Nov 2024 16:17:29 -0900
Subject: [PATCH 05/10] Adjust tests

---
 tests/test_dem.py | 15 +++------------
 xdem/dem.py       |  9 +--------
 2 files changed, 4 insertions(+), 20 deletions(-)

diff --git a/tests/test_dem.py b/tests/test_dem.py
index c6baa303..42fdcd7b 100644
--- a/tests/test_dem.py
+++ b/tests/test_dem.py
@@ -39,18 +39,11 @@ def test_init(self) -> None:
         dem3 = DEM(r)
         assert isinstance(dem3, DEM)
 
-        # From SatelliteImage
-        with warnings.catch_warnings():
-            warnings.filterwarnings("ignore", "Parse metadata from file not implemented")
-            img = gu.SatelliteImage(fn_img)
-        dem4 = DEM(img)
-        assert isinstance(dem4, DEM)
-
-        list_dem = [dem, dem2, dem3, dem4]
+        list_dem = [dem, dem2, dem3]
 
         # Check all attributes
         attrs = [at for at in _default_rio_attrs if at not in ["name", "dataset_mask", "driver"]]
-        all_attrs = attrs + gu.raster.satimg.satimg_attrs + xdem.dem.dem_attrs
+        all_attrs = attrs + xdem.dem.dem_attrs
         for attr in all_attrs:
             attrs_per_dem = [idem.__getattribute__(attr) for idem in list_dem]
             assert all(at == attrs_per_dem[0] for at in attrs_per_dem)
@@ -59,7 +52,6 @@ def test_init(self) -> None:
             (
                 np.array_equal(dem.data, dem2.data, equal_nan=True),
                 np.array_equal(dem2.data, dem3.data, equal_nan=True),
-                np.array_equal(dem3.data, dem4.data, equal_nan=True),
             )
         )
 
@@ -67,7 +59,6 @@ def test_init(self) -> None:
             (
                 np.all(dem.data.mask == dem2.data.mask),
                 np.all(dem2.data.mask == dem3.data.mask),
-                np.all(dem3.data.mask == dem4.data.mask),
             )
         )
 
@@ -192,7 +183,7 @@ def test_copy(self) -> None:
 
         # using list directly available in Class
         attrs = [at for at in _default_rio_attrs if at not in ["name", "dataset_mask", "driver"]]
-        all_attrs = attrs + gu.raster.satimg.satimg_attrs + xdem.dem.dem_attrs
+        all_attrs = attrs + xdem.dem.dem_attrs
         for attr in all_attrs:
             assert r.__getattribute__(attr) == r2.__getattribute__(attr)
 
diff --git a/xdem/dem.py b/xdem/dem.py
index 9cbe0987..168bb6e8 100644
--- a/xdem/dem.py
+++ b/xdem/dem.py
@@ -68,14 +68,7 @@ class DEM(Raster):  # type: ignore
     def __init__(
         self,
         filename_or_dataset: str | RasterType | rio.io.DatasetReader | rio.io.MemoryFile,
-        vcrs: Literal["Ellipsoid"]
-        | Literal["EGM08"]
-        | Literal["EGM96"]
-        | VerticalCRS
-        | str
-        | pathlib.Path
-        | int
-        | None = None,
+        vcrs: Literal["Ellipsoid", "EGM08", "EGM96"] | VerticalCRS | str | pathlib.Path | int | None = None,
         load_data: bool = False,
         parse_sensor_metadata: bool = False,
         silent: bool = True,

From ae3900d4e37286b9153a73f62785089a0c25b2eb Mon Sep 17 00:00:00 2001
From: Romain Hugonnet <romain.hugonnet@gmail.com>
Date: Sun, 17 Nov 2024 16:18:24 -0900
Subject: [PATCH 06/10] Last edits

---
 doc/source/conf.py | 1 -
 xdem/dem.py        | 4 ++--
 2 files changed, 2 insertions(+), 3 deletions(-)

diff --git a/doc/source/conf.py b/doc/source/conf.py
index af0066c4..7f3413ef 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -118,7 +118,6 @@
 inheritance_alias = {
     "geoutils.georaster.raster.Raster": "geoutils.Raster",
     "geoutils.georaster.raster.Mask": "geoutils.Mask",
-    "geoutils.georaster.satimg.SatelliteImage": "geoutils.SatelliteImage",
     "geoutils.geovector.Vector": "geoutils.Vector",
     "xdem.dem.DEM": "xdem.DEM",
     "xdem.coreg.base.Coreg": "xdem.Coreg",
diff --git a/xdem/dem.py b/xdem/dem.py
index 168bb6e8..d2342f24 100644
--- a/xdem/dem.py
+++ b/xdem/dem.py
@@ -81,7 +81,7 @@ def __init__(
         The vertical reference of the DEM can be defined by passing the `vcrs` argument.
         Otherwise, a vertical reference is tentatively parsed from the DEM product name.
 
-        Inherits all attributes from the :class:`geoutils.Raster` and :class:`geoutils.SatelliteImage` classes.
+        Inherits all attributes from the :class:`geoutils.Raster` class.
 
         :param filename_or_dataset: The filename of the dataset.
         :param vcrs: Vertical coordinate reference system either as a name ("WGS84", "EGM08", "EGM96"),
@@ -103,7 +103,7 @@ def __init__(
             for key in filename_or_dataset.__dict__:
                 setattr(self, key, filename_or_dataset.__dict__[key])
             return
-        # Else rely on parent SatelliteImage class options (including raised errors)
+        # Else rely on parent Raster class options (including raised errors)
         else:
             with warnings.catch_warnings():
                 warnings.filterwarnings("ignore", message="Parse metadata from file not implemented")

From 6a374695dd65dbc4a01d31c597a64ab94bf4fdb2 Mon Sep 17 00:00:00 2001
From: Romain Hugonnet <romain.hugonnet@gmail.com>
Date: Wed, 20 Nov 2024 17:39:32 -0900
Subject: [PATCH 07/10] Update geoutils version in env

---
 dev-environment.yml | 2 +-
 environment.yml     | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/dev-environment.yml b/dev-environment.yml
index 78004d13..cfaee703 100644
--- a/dev-environment.yml
+++ b/dev-environment.yml
@@ -13,7 +13,7 @@ dependencies:
   - tqdm
   - scikit-image=0.*
   - scikit-gstat>=1.0.18,<1.1
-  - geoutils=0.1.11
+  - geoutils=0.1.12
 
   # Development-specific, to mirror manually in setup.cfg [options.extras_require].
   - pip
diff --git a/environment.yml b/environment.yml
index 55e391e7..d240986f 100644
--- a/environment.yml
+++ b/environment.yml
@@ -13,7 +13,7 @@ dependencies:
   - tqdm
   - scikit-image=0.*
   - scikit-gstat>=1.0.18,<1.1
-  - geoutils=0.1.11
+  - geoutils=0.1.12
   - pip
 
    # To run CI against latest GeoUtils

From 66810266bda453dbb245293fe1e7177495b78172 Mon Sep 17 00:00:00 2001
From: Romain Hugonnet <romain.hugonnet@gmail.com>
Date: Wed, 20 Nov 2024 17:40:58 -0900
Subject: [PATCH 08/10] Proper merging

---
 requirements.txt | 2 +-
 xdem/dem.py      | 6 +-----
 2 files changed, 2 insertions(+), 6 deletions(-)

diff --git a/requirements.txt b/requirements.txt
index d8f6ecfd..fc053a11 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -11,5 +11,5 @@ scipy==1.*
 tqdm
 scikit-image==0.*
 scikit-gstat>=1.0.18,<1.1
-geoutils==0.1.11
+geoutils==0.1.12
 pip
diff --git a/xdem/dem.py b/xdem/dem.py
index 9a33feba..7b8dacf7 100644
--- a/xdem/dem.py
+++ b/xdem/dem.py
@@ -1,6 +1,3 @@
-<<<<<<< HEAD
-"""This module defines the DEM class."""
-=======
 # Copyright (c) 2024 xDEM developers
 #
 # This file is part of xDEM project:
@@ -18,8 +15,7 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-"""DEM class and functions."""
->>>>>>> upstream/main
+"""This module defines the DEM class."""
 from __future__ import annotations
 
 import pathlib

From 51c0a05a3659da9b15d5c9a5e81aac881d694af4 Mon Sep 17 00:00:00 2001
From: Romain Hugonnet <romain.hugonnet@gmail.com>
Date: Wed, 20 Nov 2024 21:15:48 -0900
Subject: [PATCH 09/10] Additional fixes

---
 doc/source/gapfill.md                                |  6 +++---
 examples/advanced/plot_demcollection.py              | 12 +++++++-----
 .../advanced/plot_variogram_estimation_modelling.py  |  1 +
 3 files changed, 11 insertions(+), 8 deletions(-)

diff --git a/doc/source/gapfill.md b/doc/source/gapfill.md
index 46bfba04..607965cb 100644
--- a/doc/source/gapfill.md
+++ b/doc/source/gapfill.md
@@ -53,9 +53,9 @@ import matplotlib.pyplot as plt
 import xdem
 
 # Load a reference DEM from 2009
-dem_2009 = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"), datetime=datetime(2009, 8, 1))
+dem_2009 = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"))
 # Load a DEM from 1990
-dem_1990 = xdem.DEM(xdem.examples.get_path("longyearbyen_tba_dem"), datetime=datetime(1990, 8, 1))
+dem_1990 = xdem.DEM(xdem.examples.get_path("longyearbyen_tba_dem"))
 # Load glacier outlines from 1990.
 glaciers_1990 = gu.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines"))
 glaciers_2010 = gu.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines_2010"))
@@ -76,7 +76,7 @@ dem_1990 = dem_1990.crop(bounds)
 We create a difference of DEMs object {class}`xdem.ddem.dDEM` to experiment on:
 
 ```{code-cell} ipython3
-ddem = xdem.dDEM(raster=dem_2009 - dem_1990, start_time=dem_1990.datetime, end_time=dem_2009.datetime)
+ddem = xdem.dDEM(raster=dem_2009 - dem_1990, start_time=datetime(1990, 8, 1), end_time=datetime(2009, 8, 1))
 
 # The example DEMs are void-free, so let's make some random voids.
 # Introduce a fifth of nans randomly throughout the dDEM.
diff --git a/examples/advanced/plot_demcollection.py b/examples/advanced/plot_demcollection.py
index c3b3fa94..58fcc750 100644
--- a/examples/advanced/plot_demcollection.py
+++ b/examples/advanced/plot_demcollection.py
@@ -2,6 +2,8 @@
 Working with a collection of DEMs
 =================================
 
+.. caution:: This functionality might be removed in future package versions.
+
 Oftentimes, more than two timestamps (DEMs) are analyzed simultaneously.
 One single dDEM only captures one interval, so multiple dDEMs have to be created.
 In addition, if multiple masking polygons exist (e.g. glacier outlines from multiple years), these should be accounted for properly.
@@ -21,8 +23,8 @@
 # We can load the DEMs as usual, but with the addition that the ``datetime`` argument should be filled.
 # Since multiple DEMs are in question, the "time dimension" is what keeps them apart.
 
-dem_2009 = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"), datetime=datetime(2009, 8, 1))
-dem_1990 = xdem.DEM(xdem.examples.get_path("longyearbyen_tba_dem"), datetime=datetime(1990, 8, 1))
+dem_2009 = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"))
+dem_1990 = xdem.DEM(xdem.examples.get_path("longyearbyen_tba_dem"))
 
 
 # %%
@@ -39,8 +41,8 @@
 
 # Fake a 2060 DEM by assuming twice the change from 1990-2009 between 2009 and 2060
 dem_2060 = dem_2009 + (dem_2009 - dem_1990).data * 3
-dem_2060.datetime = datetime(2060, 8, 1)
 
+timestamps = [datetime(1990, 8, 1), datetime(2009, 8, 1), datetime(2060, 8, 1)]
 
 # %%
 # Now, all data are ready to be collected in an :class:`xdem.DEMCollection` object.
@@ -49,8 +51,8 @@
 # 2. Two glacier outline timestamps from 1990 and 2009
 #
 
-demcollection = xdem.DEMCollection(dems=[dem_1990, dem_2009, dem_2060], outlines=outlines, reference_dem=1)
-
+demcollection = xdem.DEMCollection(dems=[dem_1990, dem_2009, dem_2060], timestamps=timestamps,
+                                   outlines=outlines, reference_dem=1)
 
 # %%
 # We can generate :class:`xdem.dDEM` objects using :func:`xdem.DEMCollection.subtract_dems`.
diff --git a/examples/advanced/plot_variogram_estimation_modelling.py b/examples/advanced/plot_variogram_estimation_modelling.py
index 5285d602..f40e49c4 100644
--- a/examples/advanced/plot_variogram_estimation_modelling.py
+++ b/examples/advanced/plot_variogram_estimation_modelling.py
@@ -39,6 +39,7 @@
 
 # %%
 # We exclude values on glacier terrain in order to isolate stable terrain, our proxy for elevation errors.
+dh.load()
 dh.set_mask(mask_glacier)
 
 # %%

From d5b01e978f742fdbfdb563e87266167f8510a680 Mon Sep 17 00:00:00 2001
From: Romain Hugonnet <romain.hugonnet@gmail.com>
Date: Wed, 20 Nov 2024 21:16:00 -0900
Subject: [PATCH 10/10] Linting

---
 examples/advanced/plot_demcollection.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/examples/advanced/plot_demcollection.py b/examples/advanced/plot_demcollection.py
index 58fcc750..8af130d3 100644
--- a/examples/advanced/plot_demcollection.py
+++ b/examples/advanced/plot_demcollection.py
@@ -51,8 +51,9 @@
 # 2. Two glacier outline timestamps from 1990 and 2009
 #
 
-demcollection = xdem.DEMCollection(dems=[dem_1990, dem_2009, dem_2060], timestamps=timestamps,
-                                   outlines=outlines, reference_dem=1)
+demcollection = xdem.DEMCollection(
+    dems=[dem_1990, dem_2009, dem_2060], timestamps=timestamps, outlines=outlines, reference_dem=1
+)
 
 # %%
 # We can generate :class:`xdem.dDEM` objects using :func:`xdem.DEMCollection.subtract_dems`.