Fix bugs following geoutils update

.github/workflows/python-package.yml

@@ -53,7 +53,7 @@ jobs:
       # pytest will run asynchronously, so test data need to be downloaded first.
     - name: Download example data
       run: |
-        python -c "from xdem.examples import *; download_longyearbyen_examples(overwrite=False)"
+        python -c "from xdem.examples import *; download_longyearbyen_examples(overwrite=True)"
     - name: Test with pytest
       run: |
         pip install pytest-cov coveralls

dev-environment.yml

@@ -12,7 +12,7 @@ dependencies:
   - matplotlib
   - opencv
   - pyproj
-  - rasterio
+  - rasterio<=1.2.10
   - scipy
   - tqdm
   - scikit-image
@@ -33,6 +33,5 @@ dependencies:
   - sphinx-gallery
 
   - pip:
-    - geoutils==0.0.8
-#    - git+https://github.com/GlacioHack/GeoUtils.git
+    - git+https://github.com/GlacioHack/geoutils.git
 

dev-requirements.txt

@@ -11,7 +11,7 @@ scipy
 rasterio
 pyproj
 tqdm
-git+https://github.com/GlacioHack/GeoUtils.git
+git+https://github.com/GlacioHack/geoutils.git
 scikit-gstat
 flake8
 opencv-contrib-python

environment.yml

@@ -18,4 +18,4 @@ dependencies:
   - proj-data
   - scikit-gstat>=0.6.8
   - pytransform3d
-  - geoutils>=0.0.5
+  - geoutils>=0.0.8

tests/test_coreg.py

@@ -199,9 +199,16 @@ def test_nuth_kaab(self):
                       transform=self.ref.transform, verbose=self.fit_params["verbose"])
 
         # Make sure that the estimated offsets are similar to what was synthesized.
-        assert abs(nuth_kaab._meta["offset_east_px"] - pixel_shift) < 0.03
-        assert abs(nuth_kaab._meta["offset_north_px"]) < 0.03
-        assert abs(nuth_kaab._meta["bias"] + bias) < 0.03
+        assert nuth_kaab._meta["offset_east_px"] == pytest.approx(pixel_shift, abs=0.03)
+        assert nuth_kaab._meta["offset_north_px"] == pytest.approx(0, abs=0.03)
+        assert nuth_kaab._meta["bias"] == pytest.approx(-bias, 0.03)
+
+        # Check that the random states forces always the same results
+        # Note: in practice, the values are not exactly equal for different OS/conda config
+        assert nuth_kaab._meta["offset_east_px"] == pytest.approx(2.000192638759735, abs=1e-7)
+        assert nuth_kaab._meta["offset_north_px"] == pytest.approx(-0.0001202906750811198, abs=1e-7)
+        assert nuth_kaab._meta["bias"] == -5.0
+
 
         # Apply the estimated shift to "revert the DEM" to its original state.
         unshifted_dem = nuth_kaab.apply(shifted_dem, transform=self.ref.transform)
@@ -566,9 +573,11 @@ def test_coreg_raster_and_ndarray_args(_) -> None:
         dem2_a = biascorr_a.apply(dem2.data, dem2.transform)
 
         # Validate that the return formats were the expected ones, and that they are equal.
+        # Issue - dem2_a does not have the same shape, the first dimension is being squeezed
+        # TODO - Fix coreg.apply?
         assert isinstance(dem2_r, xdem.DEM)
         assert isinstance(dem2_a, np.ma.masked_array)
-        assert np.array_equal(dem2_r, dem2_r)
+        assert gu.misc.array_equal(dem2_r.data.squeeze(), dem2_a, equal_nan=True)
 
         # If apply on a masked_array was given without a transform, it should fail.
         with pytest.raises(ValueError, match="'transform' must be given"):

tests/test_dem.py

@@ -4,6 +4,7 @@
 import warnings
 
 import geoutils.georaster as gr
+from geoutils.georaster.raster import _default_rio_attrs
 import geoutils.satimg as si
 import numpy as np
 import pyproj
@@ -44,7 +45,7 @@ def test_init(self):
         list_dem = [dem, dem2, dem3, dem4]
 
         # Check all attributes
-        attrs = [at for at in r._get_rio_attrs() if at not in ['name', 'dataset_mask', 'driver']]
+        attrs = [at for at in _default_rio_attrs if at not in ['name', 'dataset_mask', 'driver']]
         all_attrs = attrs + si.satimg_attrs + xdem.dem.dem_attrs
         for attr in all_attrs:
             attrs_per_dem = [idem.__getattribute__(attr) for idem in list_dem]
@@ -83,7 +84,7 @@ def test_copy(self):
         # dem_attrs = ['vref', 'vref_grid', 'ccrs']
 
         # using list directly available in Class
-        attrs = [at for at in r._get_rio_attrs() if at not in ['name', 'dataset_mask', 'driver']]
+        attrs = [at for at in _default_rio_attrs if at not in ['name', 'dataset_mask', 'driver']]
         all_attrs = attrs + si.satimg_attrs + xdem.dem.dem_attrs
         for attr in all_attrs:
             assert r.__getattribute__(attr) == r2.__getattribute__(attr)

tests/test_examples.py

@@ -0,0 +1,54 @@
+"""Functions to test the example data."""
+from __future__ import annotations
+
+import geoutils as gu
+import geoutils.spatial_tools
+from geoutils import Raster, Vector
+import numpy as np
+import pandas as pd
+import pytest
+
+import xdem
+from xdem import examples
+
+def load_examples() -> tuple[Raster, Raster, Vector, Raster]:
+    """Load example files to try coregistration methods with."""
+
+    ref_dem = Raster(examples.get_path("longyearbyen_ref_dem"))
+    tba_dem = Raster(examples.get_path("longyearbyen_tba_dem"))
+    glacier_mask = Vector(examples.get_path("longyearbyen_glacier_outlines"))
+    ddem = Raster(examples.get_path('longyearbyen_ddem'))
+
+    return ref_dem, tba_dem, glacier_mask, ddem
+
+class TestExamples:
+
+    ref_dem, tba_dem, glacier_mask, ddem = load_examples()
+
+    @pytest.mark.parametrize('rst_and_truevals',
+        [(ref_dem, np.array([868.6489, 623.42194, 180.57921, 267.30765, 601.67615], dtype=np.float32)),
+        (tba_dem, np.array([875.2358, 625.0544, 182.9936, 272.6586, 606.2897], dtype=np.float32)),
+        (ddem, np.array([-0.02423096, -0.71899414,  0.14256287,  1.1018677, -5.9209595], dtype=np.float32))])
+    def test_array_content(self, rst_and_truevals: tuple[Raster, np.ndarray]):
+        """Let's ensure the data arrays in the examples are always the same by checking randomly some values"""
+
+        rst = rst_and_truevals[0]
+        truevals = rst_and_truevals[1]
+        np.random.seed(42)
+        values = np.random.choice(rst.data.data.flatten(), size=5, replace=False)
+
+        # TODO: understand why the values differ slightly between different OS or python configurations
+        assert values == pytest.approx(truevals, abs=1e-3)
+
+    @pytest.mark.parametrize('rst_and_truenodata',
+                             [(ref_dem, 0),
+                              (tba_dem, 0),
+                              (ddem, 2316)])
+    def test_array_nodata(self, rst_and_truenodata: tuple[Raster, int]):
+        """Let's also check that the data arrays have always the same number of not finite values"""
+
+        rst = rst_and_truenodata[0]
+        truenodata = rst_and_truenodata[1]
+        mask = gu.spatial_tools.get_array_and_mask(rst)[1]
+
+        assert np.sum(mask) == truenodata

xdem/coreg.py

@@ -189,6 +189,7 @@ def residuals(parameters: tuple[float, float, float], y_values: np.ndarray, x_va
         return err
 
     # Estimate the a, b, and c parameters with least square minimisation
+    np.random.seed(seed=42)
     plsq = scipy.optimize.leastsq(func=residuals, x0=initial_guess, args=(y_medians, slice_bounds), full_output=1)
 
     a_parameter, b_parameter, c_parameter = plsq[0]
@@ -419,7 +420,8 @@ def fit(self: CoregType, reference_dem: np.ndarray | np.ma.masked_array | Raster
             transform: Optional[rio.transform.Affine] = None,
             weights: Optional[np.ndarray] = None,
             subsample: Union[float, int] = 1.0,
-            verbose: bool = False) -> CoregType:
+            verbose: bool = False,
+            random_state: None | np.random.RandomState | np.random.Generator | int = None) -> CoregType:
         """
         Estimate the coregistration transform on the given DEMs.
 
@@ -430,6 +432,7 @@ def fit(self: CoregType, reference_dem: np.ndarray | np.ma.masked_array | Raster
         :param weights: Optional. Per-pixel weights for the coregistration.
         :param subsample: Subsample the input to increase performance. <1 is parsed as a fraction. >1 is a pixel count.
         :param verbose: Print progress messages to stdout.
+        :param random_state: Random state or seed number to use for calculations (to fix random sampling during testing)
         """
 
         if weights is not None:
@@ -576,20 +579,18 @@ def apply(self, dem: np.ndarray | np.ma.masked_array | RasterType,
                 raise ValueError("Coreg method is non-rigid but has no implemented _apply_func")
 
         # Calculate final mask
-        final_mask = dem_mask + np.isnan(applied_dem)
+        final_mask = np.logical_or(dem_mask, np.isnan(applied_dem))
 
         # If the DEM was a masked_array, copy the mask to the new DEM
-        if hasattr(dem, "mask"):
-            applied_dem = np.ma.masked_array(applied_dem, mask=final_mask)  # type: ignore
-        # If the DEM was a Raster with a mask, copy the mask to the new DEM
-        elif hasattr(dem, "data") and hasattr(dem.data, "mask"):
+        if isinstance(dem, (np.ma.masked_array)):
             applied_dem = np.ma.masked_array(applied_dem, mask=final_mask)  # type: ignore
         else:
             applied_dem[final_mask] = np.nan
 
         # If the input was a Raster, return a Raster as well.
         if isinstance(dem, gu.Raster):
-            return dem.from_array(applied_dem, transform, dem.crs, nodata=dem.nodata)
+            return dem.copy(new_array=applied_dem)
+            # return dem.from_array(applied_dem, transform, dem.crs, nodata=dem.nodata)
 
         return applied_dem
 

xdem/dem.py

@@ -231,11 +231,10 @@ def to_vref(self, vref_name: str = 'EGM96', vref_grid: str | None = None):
 
         # Transform the grid
         transformer = Transformer.from_crs(ccrs_init, ccrs_dest)
-        meta = self.ds.meta
         zz = self.data
         xx, yy = self.coords(offset='center')
         zz_trans = transformer.transform(xx,yy,zz[0,:])[2]
         zz[0,:] = zz_trans
 
         # Update raster
-        self._update(metadata=meta,imgdata=zz)
+        self.data = zz

xdem/examples.py

@@ -42,6 +42,12 @@ def download_longyearbyen_examples(overwrite: bool = False):
         # print("Datasets exist")
         return
 
+    # If we ask for overwrite, also remove the processed test data
+    if overwrite:
+        for fn in list(FILEPATHS_PROCESSED.values()):
+            if os.path.exists(fn):
+                os.remove(fn)
+
     # Static commit hash to be bumped every time it needs to be.
     commit = "321f84d5a67666f45a196a31a2697e22bfaf3c59"
     # The URL from which to download the repository
@@ -93,13 +99,11 @@ def process_coregistered_examples(overwrite: bool =False):
     inlier_mask = ~glacier_mask.create_mask(reference_raster)
 
     nuth_kaab = xdem.coreg.NuthKaab()
-    nuth_kaab.fit(reference_raster.data, to_be_aligned_raster.data,
-                  inlier_mask=inlier_mask, transform=reference_raster.transform)
-    aligned_raster = nuth_kaab.apply(to_be_aligned_raster.data, transform=reference_raster.transform)
-    aligned_raster.data[~np.isfinite(aligned_raster)] = reference_raster.nodata
+    nuth_kaab.fit(reference_raster, to_be_aligned_raster,
+                  inlier_mask=inlier_mask)
+    aligned_raster = nuth_kaab.apply(to_be_aligned_raster)
 
-    diff = reference_raster - \
-           gu.Raster.from_array(aligned_raster.data, transform=reference_raster.transform, crs=reference_raster.crs, nodata=reference_raster.nodata)
+    diff = reference_raster - aligned_raster
 
     # Save it so that future calls won't need to recreate the file
     os.makedirs(os.path.dirname(FILEPATHS_PROCESSED['longyearbyen_ddem']), exist_ok=True)

xdem/spatialstats.py

@@ -914,7 +914,7 @@ def sample_empirical_variogram(values: Union[np.ndarray, RasterType], gsd: float
     # First, check all that the values provided are OK
     if isinstance(values, Raster):
         gsd = values.res[0]
-        values, mask = get_array_and_mask(values.data)
+        values, mask = get_array_and_mask(values)
     elif isinstance(values, (np.ndarray, np.ma.masked_array)):
         values, mask = get_array_and_mask(values)
     else:

xdem/terrain.py

@@ -18,29 +18,28 @@
     _has_rd = False
 
 
-def _rio_to_rda(ds: rio.DatasetReader) -> rd.rdarray:
+def _raster_to_rda(rst: RasterType) -> rd.rdarray:
     """
-    Get georeferenced richDEM array from rasterio dataset
-    :param ds: DEM
+    Get georeferenced richDEM array from geoutils.Raster
+    :param rst: DEM as raster
     :return: DEM
     """
-    arr = ds.read(1)
-    rda = rd.rdarray(arr, no_data=ds.get_nodatavals()[0])
-    rda.geotransform = ds.get_transform()
-    rda.projection = ds.get_gcps()
+    arr = rst.data.filled(rst.nodata).squeeze()
+    rda = rd.rdarray(arr, no_data=rst.nodata)
+    rda.geotransform = rst.transform.to_gdal()
 
     return rda
 
 
-def _get_terrainattr_richdem(ds: rio.DatasetReader, attribute='slope_radians') -> np.ndarray:
+def _get_terrainattr_richdem(rst: RasterType, attribute='slope_radians') -> np.ndarray:
     """
     Derive terrain attribute for DEM opened with rasterio. One of "slope_degrees", "slope_percentage", "aspect",
     "profile_curvature", "planform_curvature", "curvature" and others (see RichDEM documentation).
-    :param ds: DEM
+    :param rst: DEM as raster
     :param attribute: RichDEM terrain attribute
     :return:
     """
-    rda = _rio_to_rda(ds)
+    rda = _raster_to_rda(rst)
     terrattr = rd.TerrainAttribute(rda, attrib=attribute)
 
     return np.array(terrattr)
@@ -768,10 +767,7 @@ def get_terrain_attribute(
         for attr in attributes_using_richdem:
             attributes_requiring_surface_fit.remove(attr)
 
-        if isinstance(dem, gu.Raster):
-            # Prepare rasterio.Dataset to pass to RichDEM
-            ds = dem.ds
-        else:
+        if not isinstance(dem, gu.Raster):
             # Here, maybe we could pass the geotransform based on the resolution, and add a "default" projection as
             # this is mandated but likely not used by the rdarray format of RichDEM...
             # For now, not supported
@@ -842,7 +838,7 @@ def get_terrain_attribute(
     if make_slope:
 
         if use_richdem:
-            terrain_attributes["slope"] = _get_terrainattr_richdem(ds, attribute="slope_radians")
+            terrain_attributes["slope"] = _get_terrainattr_richdem(dem, attribute="slope_radians")
 
         else:
             if slope_method == "Horn":
@@ -863,10 +859,10 @@ def get_terrain_attribute(
 
         if use_richdem:
             # The aspect of RichDEM is returned in degrees, we convert to radians to match the others
-            terrain_attributes["aspect"] = np.deg2rad(_get_terrainattr_richdem(ds, attribute="aspect"))
+            terrain_attributes["aspect"] = np.deg2rad(_get_terrainattr_richdem(dem, attribute="aspect"))
             # For flat slopes, RichDEM returns a 90° aspect by default, while GDAL return a 180° aspect
             # We stay consistent with GDAL
-            slope_tmp = _get_terrainattr_richdem(ds, attribute="slope_radians")
+            slope_tmp = _get_terrainattr_richdem(dem, attribute="slope_radians")
             terrain_attributes["aspect"][slope_tmp == 0] = np.pi
 
         else:
@@ -914,7 +910,7 @@ def get_terrain_attribute(
     if make_curvature:
 
         if use_richdem:
-            terrain_attributes["curvature"] = _get_terrainattr_richdem(ds, attribute="curvature")
+            terrain_attributes["curvature"] = _get_terrainattr_richdem(dem, attribute="curvature")
 
         else:
             # Curvature is the second derivative of the surface fit equation.
@@ -927,7 +923,7 @@ def get_terrain_attribute(
     if make_planform_curvature:
 
         if use_richdem:
-            terrain_attributes["planform_curvature"] = _get_terrainattr_richdem(ds, attribute="planform_curvature")
+            terrain_attributes["planform_curvature"] = _get_terrainattr_richdem(dem, attribute="planform_curvature")
 
         else:
             # PLANC = 2(DH² + EG² -FGH)/(G²+H²)
@@ -953,7 +949,7 @@ def get_terrain_attribute(
     if make_profile_curvature:
 
         if use_richdem:
-            terrain_attributes["profile_curvature"] = _get_terrainattr_richdem(ds, attribute="profile_curvature")
+            terrain_attributes["profile_curvature"] = _get_terrainattr_richdem(dem, attribute="profile_curvature")
 
         else:
             # PROFC = -2(DG² + EH² + FGH)/(G²+H²)