Add dem Xarray accessor mirroring DEM class

dev-environment.yml

@@ -9,6 +9,8 @@ dependencies:
   - geopandas>=0.12.0
   - rasterio>=1.3,<2
   - scipy=1.*
+  - xarray
+  - rioxarray=0.*
   # Second-order or plus dependency on the above
   - numpy>=1,<3
   - pyproj>=3.4,<4
@@ -21,6 +23,7 @@ dependencies:
   - pip
 
   # Optional dependencies
+  - dask  # For scalable operations
   - numba=0.*  # For fast numerical operations (terrain)
   - matplotlib=3.*  # For plotting
   - scikit-learn  # For optimizations
@@ -62,4 +65,4 @@ dependencies:
       - -e ./
 
     # To run CI against latest GeoUtils
-#    - git+https://github.com/rhugonnet/geoutils.git
+#      - git+https://github.com/rhugonnet/geoutils.git

environment.yml

@@ -9,6 +9,8 @@ dependencies:
   - geopandas>=0.12.0
   - rasterio>=1.3,<2
   - scipy=1.*
+  - xarray
+  - rioxarray=0.*
   # Second-order or plus dependency on the above
   - numpy>=1,<3
   - pyproj>=3.4,<4

examples/advanced/plot_blockwise_coreg.py

@@ -22,7 +22,7 @@
 # sphinx_gallery_thumbnail_number = 2
 import matplotlib.pyplot as plt
 import numpy as np
-from geoutils.raster.distributed_computing import MultiprocConfig
+from geoutils.multiproc import MultiprocConfig
 
 import xdem
 

setup.cfg

@@ -50,6 +50,7 @@ include =
 
 [options.extras_require]
 opt =
+    dask
     numba
     matplotlib
     scikit-learn

tests/test_coreg/test_affine.py

@@ -13,7 +13,7 @@
 import rasterio as rio
 import scipy.optimize
 from geoutils import Raster, Vector
-from geoutils.raster.geotransformations import _translate
+from geoutils.raster.transformation import _translate
 from scipy.ndimage import binary_dilation
 
 from xdem import coreg, examples
@@ -57,8 +57,8 @@ class TestAffineCoreg:
     fit_args_rst_rst = dict(reference_elev=ref, to_be_aligned_elev=tba, inlier_mask=inlier_mask)
 
     # Convert DEMs to points with a bit of subsampling for speed-up
-    ref_pts = ref.to_pointcloud(data_column_name="z", subsample=50000, random_state=42).ds
-    tba_pts = ref.to_pointcloud(data_column_name="z", subsample=50000, random_state=42).ds
+    ref_pts = ref.to_pointcloud(data_column_name="z").ds
+    tba_pts = ref.to_pointcloud(data_column_name="z").ds
 
     # Raster-Point
     fit_args_rst_pts = dict(reference_elev=ref, to_be_aligned_elev=tba_pts, inlier_mask=inlier_mask)
@@ -195,11 +195,11 @@ def test_coreg_translations__synthetic(
         ref_shifted = ref.translate(shifts[0], shifts[1]) + shifts[2]
         # Convert to point cloud if input was point cloud
         if isinstance(elev_fit_args["to_be_aligned_elev"], gpd.GeoDataFrame):
-            ref_shifted = ref_shifted.to_pointcloud(data_column_name="z", subsample=50000, random_state=42).ds
+            ref_shifted = ref_shifted.to_pointcloud(data_column_name="z").ds
         elev_fit_args["to_be_aligned_elev"] = ref_shifted
 
         # Run coregistration
-        subsample_size = 50000 if coreg_method != coreg.CPD else 500
+        subsample_size = 1 if coreg_method != coreg.CPD else 500
         coreg_elev = horizontal_coreg.fit_and_apply(**elev_fit_args, subsample=subsample_size, random_state=42)
 
         # Check all fit parameters are the opposite of those used above, within a relative 1% (10% for ICP)
@@ -291,11 +291,11 @@ def test_coreg_vertical_translation__synthetic(self, fit_args: Any, vshift: floa
 
         # Convert to point cloud if input was point cloud
         if isinstance(elev_fit_args["to_be_aligned_elev"], gpd.GeoDataFrame):
-            ref_vshifted = ref_vshifted.to_pointcloud(data_column_name="z", subsample=50000, random_state=42).ds
+            ref_vshifted = ref_vshifted.to_pointcloud(data_column_name="z").ds
         elev_fit_args["to_be_aligned_elev"] = ref_vshifted
 
         # Fit the vertical shift model to the data
-        coreg_elev = vshiftcorr.fit_and_apply(**elev_fit_args, subsample=50000, random_state=42)
+        coreg_elev = vshiftcorr.fit_and_apply(**elev_fit_args)
 
         # Check that the right vertical shift was found
         assert vshiftcorr.meta["outputs"]["affine"]["shift_z"] == pytest.approx(-vshift, rel=10e-2)
@@ -390,13 +390,11 @@ def test_coreg_rigid__synthetic(
 
         # Convert to point cloud if input was point cloud
         if isinstance(elev_fit_args["to_be_aligned_elev"], gpd.GeoDataFrame):
-            ref_shifted_rotated = ref_shifted_rotated.to_pointcloud(
-                data_column_name="z", subsample=50000, random_state=42
-            ).ds
+            ref_shifted_rotated = ref_shifted_rotated.to_pointcloud(data_column_name="z").ds
         elev_fit_args["to_be_aligned_elev"] = ref_shifted_rotated
 
         # Run coregistration
-        subsample_size = 50000 if coreg_method != coreg.CPD else 500
+        subsample_size = 1 if coreg_method != coreg.CPD else 500
         coreg_elev = horizontal_coreg.fit_and_apply(**elev_fit_args, subsample=subsample_size, random_state=42)
 
         # Check that fit matrix is the invert of those used above, within a relative % for rotations
@@ -443,7 +441,7 @@ def test_coreg_rigid__synthetic(
         # Need to standardize by the elevation difference spread to avoid huge/small values close to infinity
         # Checking for 90% of variance as ICP cannot always resolve the small shifts
         # And only 10% of variance for CPD that can't resolve shifts at all
-        fac_reduc_var = 0.1 if coreg_method != coreg.CPD else 1.0
+        fac_reduc_var = 0.1 if coreg_method != coreg.CPD else 1.05
         assert np.nanvar(dh / np.nanstd(init_dh)) < fac_reduc_var
 
     @pytest.mark.parametrize(
@@ -506,7 +504,7 @@ def test_coreg_rigid__specific_args(self, rigid_coreg: coreg.Coreg) -> None:
         ref_shifted_rotated = coreg.apply_matrix(ref, matrix=matrix, centroid=centroid)
 
         # Coregister
-        subsample_size = 50000 if rigid_coreg.__class__.__name__ != "CPD" else 500
+        subsample_size = 1 if rigid_coreg.__class__.__name__ != "CPD" else 500
         rigid_coreg.fit(ref, ref_shifted_rotated, random_state=42, subsample=subsample_size)
 
     @pytest.mark.parametrize("coreg_method", [coreg.ICP, coreg.CPD, coreg.LZD])
@@ -523,7 +521,7 @@ def test_coreg_rigid__only_translation(self, coreg_method: coreg.Coreg) -> None:
         ref_shifted_rotated = coreg.apply_matrix(ref, matrix=matrix, centroid=centroid)
 
         # Run co-registration
-        subsample_size = 50000 if coreg_method != coreg.CPD else 500
+        subsample_size = 1 if coreg_method != coreg.CPD else 500
         c = coreg_method(subsample=subsample_size, only_translation=True)
         c.fit(ref, ref_shifted_rotated, random_state=42)
 
@@ -540,7 +538,7 @@ def test_coreg_rigid__only_translation(self, coreg_method: coreg.Coreg) -> None:
             assert np.allclose(invert_fit_shifts_translations[:3], shifts_rotations[:3], rtol=10e-1)
 
     @pytest.mark.parametrize("coreg_method", [coreg.ICP, coreg.CPD])
-    def test_coreg_rigid__standardize(self, coreg_method: coreg.Coreg) -> None:
+    def test_coreg_rigid__standardize(self, coreg_method: type[coreg.Coreg]) -> None:
 
         # Get reference elevation
         ref = self.ref
@@ -553,7 +551,7 @@ def test_coreg_rigid__standardize(self, coreg_method: coreg.Coreg) -> None:
         ref_shifted_rotated = coreg.apply_matrix(ref, matrix=matrix, centroid=centroid)
 
         # 1/ Run co-registration with standardization
-        subsample_size = 50000 if coreg_method != coreg.CPD else 500
+        subsample_size = 1 if coreg_method != coreg.CPD else 500
         c_std = coreg_method(subsample=subsample_size, standardize=True)
         c_std.fit(ref, ref_shifted_rotated, random_state=42)
 
@@ -610,11 +608,11 @@ def test_nuthkaab_initial_shift(self) -> None:
         # Get the coregistration method and expected shifts from the inputs
         inlier_mask = ~self.outlines.create_mask(ref)
 
-        c = coreg.NuthKaab(initial_shift=(0, 0, 0), subsample=50000)
+        c = coreg.NuthKaab(initial_shift=(0, 0, 0))
         dem_aligned_is = c.fit_and_apply(ref, tba, inlier_mask=inlier_mask, random_state=42)
         shifts_is = [c.meta["outputs"]["affine"][k] for k in ["shift_x", "shift_y", "shift_z"]]  # type: ignore
 
-        c = coreg.NuthKaab(subsample=50000)
+        c = coreg.NuthKaab()
         dem_aligned = c.fit_and_apply(ref, tba, inlier_mask=inlier_mask, random_state=42)
         shifts = [c.meta["outputs"]["affine"][k] for k in ["shift_x", "shift_y", "shift_z"]]  # type: ignore
 

tests/test_coreg/test_base.py

@@ -48,6 +48,8 @@ def assert_coreg_meta_equal(input1: Any, input2: Any) -> bool:
     """Short test function to check equality of coreg dictionary values."""
 
     # Different equality check based on input: number, callable, array, dataframe
+    if input1 is None:
+        return input2 is None
     if not isinstance(input1, type(input2)):
         return False
     elif isinstance(input1, (str, float, int, np.floating, np.integer, tuple, list)) or callable(input1):
@@ -137,7 +139,7 @@ def test_copy(self, coreg_class: Callable[[], Coreg]) -> None:
         # Make sure these don't appear in the copy
         assert corr_copy.meta != corr.meta
 
-    @pytest.mark.parametrize("subsample", [10, 10000, 0.5, 1])
+    @pytest.mark.parametrize("subsample", [10, 1000, 0.5, 1])
     def test_get_subsample_on_valid_mask(self, subsample: float | int) -> None:
         """Test the subsampling function called by all subclasses"""
 
@@ -192,8 +194,8 @@ def test_subsample(self, coreg_class: Any) -> None:
             fit_kwargs = {}
 
         # But can be overridden during fit
-        coreg_full.fit(**self.fit_params, subsample=10000, random_state=42, **fit_kwargs)
-        assert coreg_full.meta["inputs"]["random"]["subsample"] == 10000
+        coreg_full.fit(**self.fit_params, subsample=1000, random_state=42, **fit_kwargs)
+        assert coreg_full.meta["inputs"]["random"]["subsample"] == 1000
         # Check that the random state is properly set when subsampling explicitly or implicitly
         assert coreg_full.meta["inputs"]["random"]["random_state"] == 42
 
@@ -214,11 +216,11 @@ def test_subsample__pipeline(self) -> None:
         """Test that the subsample argument works as intended for pipelines"""
 
         # Check definition during instantiation
-        pipe = coreg.VerticalShift(subsample=200) + coreg.Deramp(subsample=5000)
+        pipe = coreg.VerticalShift(subsample=200) + coreg.Deramp(subsample=1000)
 
         # Check the arguments are properly defined
         assert pipe.pipeline[0].meta["inputs"]["random"]["subsample"] == 200
-        assert pipe.pipeline[1].meta["inputs"]["random"]["subsample"] == 5000
+        assert pipe.pipeline[1].meta["inputs"]["random"]["subsample"] == 1000
 
         # Check definition during fit
         pipe = coreg.VerticalShift() + coreg.Deramp()
@@ -387,12 +389,12 @@ def test_fit_and_apply(self, coreg_class: Any) -> None:
             fit_kwargs = {}
 
         # Perform fit, then apply
-        coreg_fit_then_apply.fit(**self.fit_params, subsample=10000, random_state=42, **fit_kwargs)
+        coreg_fit_then_apply.fit(**self.fit_params, **fit_kwargs)
         aligned_then = coreg_fit_then_apply.apply(elev=self.fit_params["to_be_aligned_elev"])
 
         # Perform fit and apply
         aligned_and = coreg_fit_and_apply.fit_and_apply(
-            **self.fit_params, subsample=10000, random_state=42, fit_kwargs=fit_kwargs
+            **self.fit_params, fit_kwargs=fit_kwargs
         )
 
         # Check outputs are the same: aligned raster, and metadata keys and values
@@ -415,11 +417,11 @@ def test_fit_and_apply__pipeline(self) -> None:
         coreg_fit_and_apply = coreg.NuthKaab() + coreg.Deramp()
 
         # Perform fit, then apply
-        coreg_fit_then_apply.fit(**self.fit_params, subsample=10000, random_state=42)
+        coreg_fit_then_apply.fit(**self.fit_params)
         aligned_then = coreg_fit_then_apply.apply(elev=self.fit_params["to_be_aligned_elev"])
 
         # Perform fit and apply
-        aligned_and = coreg_fit_and_apply.fit_and_apply(**self.fit_params, subsample=10000, random_state=42)
+        aligned_and = coreg_fit_and_apply.fit_and_apply(**self.fit_params)
 
         assert aligned_and.raster_equal(aligned_then, warn_failure_reason=True)
         assert list(coreg_fit_and_apply.pipeline[0].meta.keys()) == list(coreg_fit_then_apply.pipeline[0].meta.keys())
@@ -701,7 +703,7 @@ def test_pipeline(self) -> None:
 
         # Create a pipeline from two coreg methods.
         pipeline = coreg.CoregPipeline([coreg.VerticalShift(), coreg.NuthKaab()])
-        pipeline.fit(**self.fit_params, subsample=5000, random_state=42)
+        pipeline.fit(**self.fit_params)
 
         aligned_dem, _ = pipeline.apply(self.tba.data, transform=self.ref.transform, crs=self.ref.crs)
 
@@ -734,7 +736,7 @@ def test_pipeline_combinations__nobiasvar(self, coreg1: Callable[[], Coreg], cor
 
         # Create a pipeline from one affine and one biascorr methods.
         pipeline = coreg.CoregPipeline([coreg1(), coreg2()])
-        pipeline.fit(**self.fit_params, subsample=5000, random_state=42)
+        pipeline.fit(**self.fit_params)
 
         aligned_dem, _ = pipeline.apply(self.tba.data, transform=self.ref.transform, crs=self.ref.crs)
         assert aligned_dem.shape == self.ref.data.squeeze().shape
@@ -755,7 +757,7 @@ def test_pipeline_combinations__biasvar(
         # Create a pipeline from one affine and one biascorr methods
         pipeline = coreg.CoregPipeline([coreg1(), coreg.BiasCorr(**coreg2_init_kwargs)])  # type: ignore
         bias_vars = {"slope": xdem.terrain.slope(self.ref), "aspect": xdem.terrain.aspect(self.ref)}
-        pipeline.fit(**self.fit_params, bias_vars=bias_vars, subsample=5000, random_state=42)
+        pipeline.fit(**self.fit_params, bias_vars=bias_vars)
 
         aligned_dem, _ = pipeline.apply(
             self.tba.data, transform=self.ref.transform, crs=self.ref.crs, bias_vars=bias_vars
@@ -810,7 +812,7 @@ def test_pipeline__errors(self) -> None:
     def test_pipeline_pts(self) -> None:
 
         pipeline = coreg.NuthKaab() + coreg.DhMinimize()
-        ref_points = self.ref.to_pointcloud(subsample=5000, random_state=42)
+        ref_points = self.ref.to_pointcloud()
 
         # Check that this runs without error
         pipeline.fit(reference_elev=ref_points, to_be_aligned_elev=self.tba)

tests/test_coreg/test_biascorr.py

@@ -41,8 +41,8 @@ class TestBiasCorr:
     fit_args_rst_rst = dict(reference_elev=ref, to_be_aligned_elev=tba, inlier_mask=inlier_mask)
 
     # Convert DEMs to points with a bit of subsampling for speed-up
-    tba_pts = tba.to_pointcloud(data_column_name="z", subsample=50000, random_state=42)
-    ref_pts = ref.to_pointcloud(data_column_name="z", subsample=50000, random_state=42)
+    tba_pts = tba.to_pointcloud(data_column_name="z")
+    ref_pts = ref.to_pointcloud(data_column_name="z")
 
     # Raster-Point
     fit_args_rst_pts = dict(reference_elev=ref, to_be_aligned_elev=tba_pts, inlier_mask=inlier_mask)
@@ -292,7 +292,7 @@ def test_biascorr__bin_2d(self, fit_args: Any, bin_sizes: Any, bin_statistic: An
         elev_fit_args.update({"bias_vars": bias_vars_dict})
 
         # Run with input parameter, and using only 100 subsamples for speed
-        bcorr.fit(**elev_fit_args, subsample=10000, random_state=42)
+        bcorr.fit(**elev_fit_args)
 
         # Check that variable names are defined during fit
         assert bcorr.meta["inputs"]["fitorbin"]["bias_var_names"] == ["elevation", "slope"]
@@ -439,7 +439,7 @@ def test_directionalbias__synthetic(self, fit_args: Any, angle: float, nb_freq:
             plt.show()
 
             dirbias = biascorr.DirectionalBias(angle=angle, fit_or_bin="bin", bin_sizes=10000)
-            dirbias.fit(reference_elev=self.ref, to_be_aligned_elev=bias_dem, subsample=10000, random_state=42)
+            dirbias.fit(reference_elev=self.ref, to_be_aligned_elev=bias_dem)
             xdem.spatialstats.plot_1d_binning(
                 df=dirbias.meta["outputs"]["fitorbin"]["bin_dataframe"],
                 var_name="angle",
@@ -464,14 +464,12 @@ def test_directionalbias__synthetic(self, fit_args: Any, angle: float, nb_freq:
         elev_fit_args = fit_args.copy()
         if isinstance(elev_fit_args["to_be_aligned_elev"], gpd.GeoDataFrame):
             # Need a higher sample size to get the coefficients right here
-            bias_elev = bias_dem.to_pointcloud(data_column_name="z", subsample=50000, random_state=42).ds
+            bias_elev = bias_dem.to_pointcloud(data_column_name="z").ds
         else:
             bias_elev = bias_dem
         dirbias.fit(
             elev_fit_args["reference_elev"],
             to_be_aligned_elev=bias_elev,
-            subsample=40000,
-            random_state=42,
             bounds_amp_wave_phase=bounds,
             niter=2,
         )
@@ -524,10 +522,10 @@ def test_deramp__synthetic(self, fit_args: Any, order: int) -> None:
         deramp = biascorr.Deramp(poly_order=order)
         elev_fit_args = fit_args.copy()
         if isinstance(elev_fit_args["to_be_aligned_elev"], gpd.GeoDataFrame):
-            bias_elev = bias_dem.to_pointcloud(data_column_name="z", subsample=30000, random_state=42).ds
+            bias_elev = bias_dem.to_pointcloud(data_column_name="z").ds
         else:
             bias_elev = bias_dem
-        deramp.fit(elev_fit_args["reference_elev"], to_be_aligned_elev=bias_elev, subsample=20000, random_state=42)
+        deramp.fit(elev_fit_args["reference_elev"], to_be_aligned_elev=bias_elev)
 
         # Check high-order fit parameters are the same within 10%
         fit_params = deramp.meta["outputs"]["fitorbin"]["fit_params"]
@@ -582,14 +580,12 @@ def test_terrainbias__synthetic(self, fit_args: Any) -> None:
         )
         elev_fit_args = fit_args.copy()
         if isinstance(elev_fit_args["to_be_aligned_elev"], gpd.GeoDataFrame):
-            bias_elev = bias_dem.to_pointcloud(data_column_name="z", subsample=20000, random_state=42).ds
+            bias_elev = bias_dem.to_pointcloud(data_column_name="z").ds
         else:
             bias_elev = bias_dem
         tb.fit(
             elev_fit_args["reference_elev"],
             to_be_aligned_elev=bias_elev,
-            subsample=10000,
-            random_state=42,
             bias_vars={"max_curvature": maxc},
         )
 

tests/test_coreg/test_blockwise.py

@@ -11,8 +11,7 @@
 import pytest
 from geoutils import Raster, Vector
 from geoutils.interface.gridding import _grid_pointcloud
-from geoutils.raster import ClusterGenerator
-from geoutils.raster.distributed_computing import MultiprocConfig
+from geoutils.multiproc import MultiprocConfig, ClusterGenerator
 
 import xdem
 from xdem.coreg import BlockwiseCoreg, Coreg