Update with SatelliteImage class removal and restructuration in GeoUtils

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.12
 
   # Development-specific, to mirror manually in setup.cfg [options.extras_require].
   - pip

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",

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.

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

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.12
   - pip
 
    # To run CI against latest GeoUtils

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,9 @@
 # 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`.

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)
 
 # %%

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.12
 pip

tests/test_coreg/test_affine.py

@@ -13,7 +13,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
@@ -129,7 +129,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
         )

tests/test_dem.py

@@ -40,18 +40,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)
@@ -60,15 +53,13 @@ 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),
             )
         )
 
         assert np.logical_and.reduce(
             (
                 np.all(dem.data.mask == dem2.data.mask),
                 np.all(dem2.data.mask == dem3.data.mask),
-                np.all(dem3.data.mask == dem4.data.mask),
             )
         )
 
@@ -193,7 +184,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)
 

xdem/coreg/affine.py

@@ -35,8 +35,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

xdem/coreg/base.py

@@ -49,17 +49,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
@@ -540,7 +540,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):
@@ -1263,7 +1263,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
@@ -1279,7 +1279,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

xdem/ddem.py

@@ -26,7 +26,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
 

xdem/dem.py

@@ -15,7 +15,7 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-"""DEM class and functions."""
+"""This module defines the DEM class."""
 from __future__ import annotations
 
 import pathlib
@@ -26,7 +26,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
@@ -52,7 +52,7 @@
 dem_attrs = ["_vcrs", "_vcrs_name", "_vcrs_grid"]
 
 
-class DEM(SatelliteImage):  # type: ignore
+class DEM(Raster):  # type: ignore
     """
     The digital elevation model.
 
@@ -85,24 +85,29 @@ class DEM(SatelliteImage):  # 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,
-        **kwargs: Any,
+        downsample: int = 1,
+        nodata: int | float | None = None,
     ) -> None:
         """
         Instantiate a digital elevation model.
 
         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"),
             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
@@ -115,11 +120,18 @@ 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")
-                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:
@@ -146,8 +158,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:
@@ -199,7 +211,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,

xdem/spatialstats.py

@@ -33,14 +33,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

xdem/volume.py

@@ -27,8 +27,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,