Merge a5cffee into 057b8c4

.gitignore

@@ -153,6 +153,7 @@ examples/data/*.csv
 doc/source/basic_examples/
 doc/source/advanced_examples/
 doc/source/gen_modules/
+doc/source/sg_execution_times.rst
 examples/basic/temp.tif
 
 # Directory where myst_nb executes jupyter code

dev-environment.yml

@@ -3,7 +3,7 @@ channels:
   - conda-forge
 dependencies:
   - python>=3.9,<3.12
-  - geopandas>=0.12.0,<0.14
+  - geopandas>=0.12.0
   - numba=0.*
   - numpy=1.*
   - matplotlib=3.*
@@ -13,7 +13,7 @@ dependencies:
   - tqdm
   - scikit-image=0.*
   - scikit-gstat>=1.0
-  - geoutils=0.0.15
+#  - geoutils=0.0.15
 
   # Development-specific, to mirror manually in setup.cfg [options.extras_require].
   - pip
@@ -51,4 +51,4 @@ dependencies:
     - noisyopt
 
     # To run CI against latest GeoUtils
-#    - git+https://github.com/GlacioHack/geoutils.git
+    - git+https://github.com/GlacioHack/geoutils.git

doc/source/about_xdem.md

@@ -20,7 +20,7 @@ Additionally, xDEM aims to be **efficient**, **scalable** and **state-of-the-art
 :class: margin
 xDEM is in early stages of development and its features might evolve rapidly. Note the version you are working on for
 **reproducibility**!
-We are working on making features fully consistent for the first long-term release ``v0.1`` (likely sometime in 2023).
+We are working on making features fully consistent for the first long-term release ``v0.1`` (planned early 2024).
 ```
 
 In details, those mean:

doc/source/index.md

@@ -74,7 +74,7 @@ Dive into the full documentation.
 :::{important}
 xDEM is in early stages of development and its features might evolve rapidly. Note the version you are
 working on for reproducibility!
-We are working on making features fully consistent for the first long-term release `v0.1` (likely sometime in 2023).
+We are working on making features fully consistent for the first long-term release `v0.1` (planned early 2024).
 :::
 
 ```{toctree}

doc/source/quick_start.md

@@ -1,49 +1,104 @@
+---
+file_format: mystnb
+jupytext:
+  formats: md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+kernelspec:
+  display_name: xdem-env
+  language: python
+  name: xdem
+---
 (quick-start)=
 
 # Quick start
 
-## Sample data
+Below is a short example show-casing some of the core functionalities of xDEM.
+To find an example about a specific functionality, jump directly to {ref}`quick-gallery`.
 
-xDEM comes with some sample data that is used throughout this documentation to demonstrate the features. If not done already, the sample data can be downloaded with the command
+## Short example
 
-```python
-xdem.examples.download_longyearbyen_examples()
+```{note}
+:class: margin
+xDEM relies largely on [its sister-package GeoUtils](https://geoutils.readthedocs.io/) for geospatial handling 
+(reprojection, cropping, raster-vector interface, point interpolation) as well as numerics 
+(NumPy interface). 🙂
 ```
 
-The dataset keys and paths can be found from
+xDEM revolves around the {class}`~xdem.DEM` class (a subclass of {class}`~geoutils.Raster`), from 
+which most methods can be called and the {class}`~xdem.coreg.Coreg` classes to build modular coregistration pipelines.
 
-```python
-xdem.examples.available
-```
-
-## Load DEMs and calculate elevation difference
+Below, in a few lines, we load two DEMs and a vector of glacier outlines, crop them to a common extent, 
+align the DEMs using coregistration and estimate a map of error in elevation difference using stable terrain, and 
+finally plot and save the result!
 
-```python
+```{code-cell} ipython3
 import xdem
+import geoutils as gu
+
+# Examples files: paths to two GeoTIFFs and an ESRI shapefile of glacier outlines
+filename_dem_ref = xdem.examples.get_path("longyearbyen_ref_dem")
+filename_dem_tba = xdem.examples.get_path("longyearbyen_tba_dem")
+filename_glacier_outlines = xdem.examples.get_path("longyearbyen_glacier_outlines")
+
+# Open files by instantiating DEM and Vector
+# (DEMs are loaded lazily = only metadata but not array unless required)
+dem_ref = xdem.DEM(filename_dem_ref)
+dem_tba = xdem.DEM(filename_dem_tba)
+vect_gla = gu.Vector(filename_glacier_outlines)
+
+# Clip outlines to extent of reference DEM (method from GeoUtils)
+vect_gla = vect_gla.crop(dem_ref, clip=True)
+
+# Create a mask from glacier polygons (method from GeoUtils)
+mask_gla = vect_gla.create_mask(dem_ref)
+
+# We convert the vertical CRS of one DEM to the EGM96 geoid
+dem_ref.to_vcrs("EGM96", force_source_vcrs="Ellipsoid")
+
+# Align the two DEMs with a coregistration: 3D shift + 2nd-order 2D poly
+mycoreg = xdem.coreg.NuthKaab() + xdem.coreg.Deramp(poly_order=2)
+mycoreg.fit(dem_ref, dem_tba, inlier_mask=~mask_gla)
+dem_aligned = mycoreg.apply(dem_tba)
 
-# Load data
-dem_2009 = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"))
-dem_1990 = xdem.DEM(xdem.examples.get_path("longyearbyen_tba_dem_coreg"))
+# Get elevation difference
+dh = dem_ref - dem_aligned
 
-# Calculate the difference
-ddem = dem_2009 - dem_1990
+# Derive slope and curvature attributes
+slope, maximum_curvature = xdem.terrain.get_terrain_attribute(
+    dem_ref, attribute=["slope", "maximum_curvature"]
+)
 
-# Plot
-ddem.show(cmap='coolwarm_r', vmin=-20, vmax=20, cb_title="Elevation change (m)")
+# Estimate elevation change error from stable terrain as a function of slope and curvature 
+dem_error = xdem.spatialstats.infer_heteroscedasticity_from_stable(
+    dh, list_var=[slope, maximum_curvature], unstable_mask=mask_gla
+)[0]
+
+# Plot error map of elevation difference and glacier outlines
+dem_error.show(vmin=2, vmax=7, cmap="Reds", cbar_title=r"Elevation error (1$\sigma$, m)")
+vect_gla.show(dem_error, fc='none', ec='k', lw=0.5)
 
 # Save to file
-ddem.save("temp.tif")
+dem_error.save("my_dem_error.tif")
+```
+
+```{code-cell} ipython3
+:tags: [remove-cell]
+import os
+os.remove("my_dem_error.tif")
 ```
 
-A detailed example on how to load raster data, reproject it to a different grid/projection, run simple arithmetic operations such as subtraction, plotting the data and saving to file can be found in the example gallery {ref}`sphx_glr_basic_examples_plot_dem_subtraction.py`.
-
-% .. raw:: html
-%
-%     <div class="sphx-glr-thumbcontainer" tooltip="DEM subtraction">
-%
-% .. only:: html
-%
-%  .. figure:: /auto_examples/images/thumb/sphx_glr_plot_dem_subtraction_thumb.png
-%      :alt: DEM subtraction
-%
-%      :ref:`sphx_glr_auto_examples_plot_dem_subtraction.py`
+(quick-gallery)=
+## More examples
+
+To dive into more illustrated code, explore our gallery of examples that is composed of:
+- An {ref}`examples-basic` section on simpler routines (terrain attributes, pre-defined coregistration and uncertainty pipelines),
+- An {ref}`examples-advanced` section using advanced pipelines (for in-depth coregistration and uncertainty analysis).
+
+See also the concatenated list of examples below.
+
+```{eval-rst}
+.. minigallery:: xdem.DEM
+    :add-heading: Examples using DEMs
+```

doc/source/spatialstats.md

@@ -175,7 +175,7 @@ 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.raster.subsample_array(dh_arr, subsample=10000, return_indices=True, random_state=42)
+indices = gu.raster.sample_array(dh_arr, sample=10000, return_indices=True, random_state=42)
 dh_arr = dh_arr[indices]
 slope_arr = slope_arr[indices]
 ```

environment.yml

@@ -3,7 +3,7 @@ channels:
   - conda-forge
 dependencies:
   - python>=3.9,<3.12
-  - geopandas>=0.12.0,<0.14
+  - geopandas>=0.12.0
   - numba=0.*
   - numpy=1.*
   - matplotlib=3.*
@@ -13,8 +13,8 @@ dependencies:
   - tqdm
   - scikit-image=0.*
   - scikit-gstat>=1.0
-  - geoutils=0.0.15
+#  - geoutils=0.0.15
   - pip
 
-#  - pip:
-#     - git+https://github.com/GlacioHack/geoutils.git
+  - pip:
+     - git+https://github.com/GlacioHack/geoutils.git

examples/advanced/README.rst

@@ -1,2 +1,7 @@
+.. _examples-advanced:
+
 Advanced
 ========
+
+Examples for setting up **specific coregistration or bias-correction pipelines**, **comparing terrain methods**,
+or **refining an error model for DEM uncertainty analysis**.

examples/basic/README.rst

@@ -1,2 +1,7 @@
+.. _examples-basic:
+
 Basic
 =====
+
+Exmaples using **terrain methods** and **DEM differences**, as well as
+pre-defined **coregistration** and **uncertainty analysis** pipelines.

examples/basic/plot_icp_coregistration.py

@@ -22,7 +22,7 @@
 dem = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"))
 
 subset_extent = [523000, 8660000, 529000, 8665000]
-dem.crop(subset_extent)
+dem = dem.crop(subset_extent)
 
 # %%
 # Let's plot a hillshade of the mountain for context.

requirements.txt

@@ -1,7 +1,7 @@
 # This file is auto-generated from environment.yml, do not modify.
 # See that file for comments about the need/usage of each dependency.
 
-geopandas>=0.12.0,<0.14
+geopandas>=0.12.0
 numba==0.*
 numpy==1.*
 matplotlib==3.*
@@ -11,7 +11,7 @@ scipy==1.*
 tqdm
 scikit-image==0.*
 scikit-gstat>=1.0
-geoutils==0.0.15
 pip
+geoutils
 setuptools>=64
 setuptools_scm[toml]>=8

tests/test_coreg/test_affine.py

@@ -182,9 +182,9 @@ def test_coreg_example_shift(self, shift_px, coreg_class, points_or_raster, verb
 
         # shift DEM by shift_px
         shifted_ref = self.ref.copy()
-        shifted_ref.shift(shift_px[0] * res, shift_px[1] * res)
+        shifted_ref.shift(shift_px[0] * res, shift_px[1] * res, inplace=True)
 
-        shifted_ref_points = shifted_ref.to_points(as_array=False, subset=subsample, pixel_offset="center").ds
+        shifted_ref_points = shifted_ref.to_points(as_array=False, sample=subsample, pixel_offset="center").ds
         shifted_ref_points["E"] = shifted_ref_points.geometry.x
         shifted_ref_points["N"] = shifted_ref_points.geometry.y
         shifted_ref_points.rename(columns={"b1": "z"}, inplace=True)

tests/test_coreg/test_base.py

@@ -259,7 +259,7 @@ def test_coreg_raster_and_ndarray_args(self) -> None:
         dem1.data[1, 1] = 100
 
         # Translate the DEM 1 "meter" right and add a vertical shift
-        dem2 = dem1.reproject(dst_bounds=rio.coords.BoundingBox(1, 0, 6, 5), silent=True)
+        dem2 = dem1.reproject(bounds=rio.coords.BoundingBox(1, 0, 6, 5), silent=True)
         dem2 += 1
 
         # Create a vertical shift correction for Rasters ("_r") and for arrays ("_a")
@@ -374,7 +374,7 @@ def test_apply_resample(self, inputs: list[Any]) -> None:
                 "dem1.crs",
                 "fit",
                 "warns",
-                "'reference_dem' .* overrides the given 'transform'",
+                "'reference_dem' .* overrides the given *",
             ),
             ("dem1.data", "dem2", "dem1.transform", "None", "fit", "warns", "'dem_to_be_aligned' .* overrides .*"),
             (
@@ -635,7 +635,7 @@ def test_pipeline_pts(self) -> None:
         warnings.simplefilter("ignore")
 
         pipeline = coreg.NuthKaab() + coreg.GradientDescending()
-        ref_points = self.ref.to_points(as_array=False, subset=5000, pixel_offset="center").ds
+        ref_points = self.ref.to_points(as_array=False, sample=5000, pixel_offset="center").ds
         ref_points["E"] = ref_points.geometry.x
         ref_points["N"] = ref_points.geometry.y
         ref_points.rename(columns={"b1": "z"}, inplace=True)
@@ -783,8 +783,8 @@ def test_blockwise_coreg(self, pipeline: Coreg, subdivision: int) -> None:
     def test_blockwise_coreg_large_gaps(self) -> None:
         """Test BlockwiseCoreg when large gaps are encountered, e.g. around the frame of a rotated DEM."""
         warnings.simplefilter("error")
-        reference_dem = self.ref.reproject(dst_crs="EPSG:3413", dst_res=self.ref.res, resampling="bilinear")
-        dem_to_be_aligned = self.tba.reproject(dst_ref=reference_dem, resampling="bilinear")
+        reference_dem = self.ref.reproject(crs="EPSG:3413", res=self.ref.res, resampling="bilinear")
+        dem_to_be_aligned = self.tba.reproject(ref=reference_dem, resampling="bilinear")
 
         blockwise = xdem.coreg.BlockwiseCoreg(xdem.coreg.NuthKaab(), 64, warn_failures=False)
 

tests/test_dem.py

@@ -98,7 +98,7 @@ def test_init__vcrs(self) -> None:
         # Tests 2: instantiation with a file that has a 3D CRS
         # Create such a file
         dem = DEM(fn_img)
-        dem_reproj = dem.reproject(dst_crs=4979)
+        dem_reproj = dem.reproject(crs=4979)
 
         # Save to temporary folder
         temp_dir = tempfile.TemporaryDirectory()
@@ -186,7 +186,7 @@ def test_copy(self) -> None:
         assert isinstance(r2, xdem.dem.DEM)
 
         # Check all immutable attributes are equal
-        # raster_attrs = ['bounds', 'count', 'crs', 'dtypes', 'height', 'indexes', 'nodata',
+        # raster_attrs = ['bounds', 'count', 'crs', 'dtypes', 'height', 'bands', 'nodata',
         #                    'res', 'shape', 'transform', 'width']
         # satimg_attrs = ['satellite', 'sensor', 'product', 'version', 'tile_name', 'datetime']
         # dem_attrs = ['vcrs', 'vcrs_grid', 'vcrs_name', 'ccrs']
@@ -263,15 +263,15 @@ def test_to_vcrs(self) -> None:
         dem = DEM(fn_dem)
 
         # Reproject in WGS84 2D
-        dem = dem.reproject(dst_crs=4326)
+        dem = dem.reproject(crs=4326)
         dem_before_trans = dem.copy()
 
         # Set ellipsoid as vertical reference
         dem.set_vcrs(new_vcrs="Ellipsoid")
         ccrs_init = dem.ccrs
         median_before = np.nanmean(dem)
         # Transform to EGM96 geoid
-        dem.to_vcrs(dst_vcrs="EGM96")
+        dem.to_vcrs(vcrs="EGM96")
         median_after = np.nanmean(dem)
 
         # About 32 meters of difference in Svalbard between EGM96 geoid and ellipsoid

tests/test_spatialstats.py

@@ -72,9 +72,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.raster.subsample_array(
-            self.diff.data.flatten(), subsample=10000, return_indices=True, random_state=42
-        )
+        indices = gu.raster.sample_array(self.diff.data.flatten(), sample=10000, return_indices=True, random_state=42)
 
         # 1D binning, by default will create 10 bins
         df = xdem.spatialstats.nd_binning(

xdem/coreg/base.py

@@ -42,8 +42,8 @@
     RasterType,
     get_array_and_mask,
     raster,
+    sample_array,
     subdivide_array,
-    subsample_array,
 )
 from tqdm import tqdm
 
@@ -450,7 +450,7 @@ def residuals(coefs: NDArrayf, x_coords: NDArrayf, y_coords: NDArrayf, targets:
         print("Estimating deramp function...")
 
     # reduce number of elements for speed
-    rand_indices = subsample_array(x_coords, subsample=subsample, return_indices=True)
+    rand_indices = sample_array(x_coords, sample=subsample, return_indices=True)
     x_coords = x_coords[rand_indices]
     y_coords = y_coords[rand_indices]
     ddem = ddem[rand_indices]
@@ -748,9 +748,9 @@ def _get_subsample_on_valid_mask(self, valid_mask: NDArrayb, verbose: bool = Fal
             # Build a low memory masked array with invalid values masked to pass to subsampling
             ma_valid = np.ma.masked_array(data=np.ones(np.shape(valid_mask), dtype=bool), mask=~valid_mask)
             # Take a subsample within the valid values
-            indices = gu.raster.subsample_array(
+            indices = gu.raster.sample_array(
                 ma_valid,
-                subsample=self._meta["subsample"],
+                sample=self._meta["subsample"],
                 return_indices=True,
                 random_state=self._meta["random_state"],
             )
@@ -1041,8 +1041,8 @@ def fit_pts(
         if subsample != 1.0:
 
             # Randomly pick N inliers in the full_mask where N=subsample
-            random_valids = subsample_array(
-                ref_dem[z_name].values, subsample=subsample, return_indices=True, random_state=random_state
+            random_valids = sample_array(
+                ref_dem[z_name].values, sample=subsample, return_indices=True, random_state=random_state
             )
 
             # Subset to the N random inliers