From 54657dff55c8f52612a5075325d0dc1a98dada10 Mon Sep 17 00:00:00 2001
From: adehecq <3285905+adehecq@users.noreply.github.com>
Date: Thu, 12 May 2022 17:04:36 +0200
Subject: [PATCH 01/11] Remove unused imports

---
 xdem/coreg.py | 5 -----
 1 file changed, 5 deletions(-)

diff --git a/xdem/coreg.py b/xdem/coreg.py
index a6a754c8..cbf30293 100644
--- a/xdem/coreg.py
+++ b/xdem/coreg.py
@@ -3,12 +3,7 @@
 
 import copy
 import concurrent.futures
-import json
-import os
-import subprocess
-import tempfile
 import warnings
-from enum import Enum
 from typing import Any, Callable, Optional, overload, Union, Sequence, TypeVar
 
 try:

From 25271a759c97601716dcd0f946a17d1006c7f11c Mon Sep 17 00:00:00 2001
From: adehecq <3285905+adehecq@users.noreply.github.com>
Date: Thu, 12 May 2022 17:06:33 +0200
Subject: [PATCH 02/11] Include a one-liner function for DEM coregistration

---
 xdem/coreg.py | 119 ++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 119 insertions(+)

diff --git a/xdem/coreg.py b/xdem/coreg.py
index cbf30293..5870f4d4 100644
--- a/xdem/coreg.py
+++ b/xdem/coreg.py
@@ -15,6 +15,8 @@
 import geoutils as gu
 from geoutils.georaster import RasterType
 from geoutils import spatial_tools
+from geoutils._typing import AnyNumber
+import matplotlib.pyplot as plt
 import numpy as np
 import rasterio as rio
 import rasterio.warp  # pylint: disable=unused-import
@@ -1979,3 +1981,120 @@ def warp_dem(
     assert not np.all(np.isnan(warped)), "All-NaN output."
 
     return warped.reshape(dem.shape)
+
+
+def dem_coregistration(
+    src_dem_path: str,
+    ref_dem_path: str,
+    shpfile: str,
+    out_dem_path: str,
+    coreg_method: Coreg | None = NuthKaab() + BiasCorr(bias_func=np.nanmedian),
+    grid: str = "ref",
+    filtering: bool = True,
+    slope_lim: list(AnyNumber) = (0.1, 40),
+    plot: bool = False,
+    out_fig: str = None,
+    verbose: bool = False,
+):
+    """
+    Coregister a selected DEM to a reference DEM.
+    Reads both DEMs, reproject DEM onto ref DEM grid, mask content of shpfile, run the coregistration and save the coregistered DEM as well as some optional figures and returns some statistics.
+
+    :param src_dem_path: path to the input DEM to be coregistered
+    :param ref_dem: path to the reference DEM
+    :param shpfile: path to a vector file containing areas to be masked for coregistration
+    :param out_dem_path: Path where to save the coregistered DEM
+    :param coreg_method: The xdem coregistration method, or pipeline. If set to None, DEMs will be resampled to ref grid and optionally filtered, but not coregistered.
+    :param grid: the grid to be used during coregistration, set either to "ref" or "src".
+    :param filtering: if set to True, filtering will be applied prior to coregistration
+    :param plot: Set to True to plot a figure of elevation diff before/after coregistration
+    :param out_fig: Path to the output figure. If None will display to screen.
+    :param verbose: set to True to print details on screen during coregistration.
+
+    :returns: a tuple containing - basename of coregistered DEM, [count of obs, median and NMAD over stable terrain, coverage over roi] before coreg, [same stats] after coreg
+    """
+    # Load both DEMs
+    ref_dem = xdem.DEM(ref_dem_path)
+    src_dem = xdem.DEM(src_dem_path)
+
+    # Reproject to common grid
+    if grid == "ref":
+        src_dem = src_dem.reproject(ref_dem, resampling='bilinear', silent=True)
+    elif grid == "src":
+        ref_dem = ref_dem.reproject(src_dem, resampling='bilinear', silent=True)
+    else:
+        raise ValueError(f"`grid` must be either 'ref' or 'src' - currently set to {grid}")
+
+    # Create raster mask
+    outlines = gu.Vector(shpfile)
+    stable_mask = ~outlines.create_mask(src_dem)
+
+    # Calculate dDEM
+    ddem = src_dem - ref_dem
+
+    # Filter gross outliers in stable terrain
+    if filtering:
+        # TO DO implement the NMAD filter in xdem
+        inlier_mask = stable_mask  # nmad_filter(ddem.data, stable_mask, verbose=False)
+
+        # Exclude steep slopes for coreg
+        slope = xdem.terrain.slope(ref_dem)
+        inlier_mask[slope.data < slope_lim[0]] = False
+        inlier_mask[slope.data > slope_lim[1]] = False
+
+    else:
+        inlier_mask = stable_mask
+
+    # Calculate dDEM statistics on pixels used for coreg
+    inlier_data = ddem.data[inlier_mask].compressed()
+    nstable_orig, med_orig, nmad_orig = len(inlier_data), np.median(inlier_data), xdem.spatialstats.nmad(inlier_data)
+
+    # Coregister to reference - Note: this will spread NaN
+    # Better strategy: calculate shift, update transform, resample
+    if isinstance(coreg_method, xdem.coreg.Coreg):
+        coreg_method.fit(ref_dem, src_dem, inlier_mask, verbose=verbose)
+        dem_coreg = coreg_method.apply(src_dem, dilate_mask=False)
+    elif coreg_method is None:
+        dem_coreg = src_dem
+    ddem_coreg = dem_coreg - ref_dem
+
+    # Calculate new stats
+    inlier_data = ddem_coreg.data[inlier_mask].compressed()
+    nstable_coreg, med_coreg, nmad_coreg = len(inlier_data), np.median(inlier_data), xdem.spatialstats.nmad(inlier_data)
+
+    # Plot results
+    if plot:
+        # Max colorbar value - 98th percentile rounded to nearest 5
+        vmax = np.percentile(np.abs(ddem.data.compressed()), 98) // 5 * 5
+
+        plt.figure(figsize=(11, 5))
+
+        ax1 = plt.subplot(121)
+        plt.imshow(ddem.data.squeeze(), cmap="coolwarm_r", vmin=-vmax, vmax=vmax)
+        cb = plt.colorbar()
+        cb.set_label("Elevation change (m)")
+        ax1.set_title(f"Before coreg\n\nmed = {med_orig:.2f} m - NMAD = {nmad_orig:.2f} m")
+
+        ax2 = plt.subplot(122, sharex=ax1, sharey=ax1)
+        plt.imshow(ddem_coreg.data.squeeze(), cmap="coolwarm_r", vmin=-vmax, vmax=vmax)
+        cb = plt.colorbar()
+        cb.set_label("Elevation change (m)")
+        ax2.set_title(f"After coreg\n\nmed = {med_coreg:.2f} m - NMAD = {nmad_coreg:.2f} m")
+
+        plt.tight_layout()
+        if out_fig is None:
+            plt.show()
+        else:
+            plt.savefig(out_fig, dpi=200)
+            plt.close()
+
+    # Save coregistered DEM
+    dem_coreg.save(out_dem_path, tiled=True)
+
+    # Save stats to DataFrame
+    out_stats = pd.DataFrame(
+        ((nstable_orig, med_orig, nmad_orig, nstable_coreg, med_coreg, nmad_coreg),),
+        columns=("nstable_orig", "med_orig", "nmad_orig", "nstable_coreg", "med_coreg", "nmad_coreg")
+    )
+
+    return dem_coreg, out_stats

From 81430bf10e3071a183da6b83047d7be84a71e674 Mon Sep 17 00:00:00 2001
From: adehecq <3285905+adehecq@users.noreply.github.com>
Date: Fri, 20 May 2022 16:31:06 +0200
Subject: [PATCH 03/11] Split coreg into horizontal/vertical coreg. Update
 stats and plots

---
 xdem/coreg.py | 66 +++++++++++++++++++++++++++++++++++++++++++--------
 1 file changed, 56 insertions(+), 10 deletions(-)

diff --git a/xdem/coreg.py b/xdem/coreg.py
index 5870f4d4..e8c854b8 100644
--- a/xdem/coreg.py
+++ b/xdem/coreg.py
@@ -1983,12 +1983,27 @@ def warp_dem(
     return warped.reshape(dem.shape)
 
 
+hmodes_dict = {
+    "nuth_kaab": NuthKaab(),
+    "nuth_kaab_block": BlockwiseCoreg(coreg=NuthKaab(), subdivision=16),
+    "icp": ICP(),
+}
+
+vmodes_dict = {
+    "median": BiasCorr(bias_func=np.median),
+    "mean": BiasCorr(bias_func=np.mean),
+    "deramp": Deramp(),
+}
+
 def dem_coregistration(
     src_dem_path: str,
     ref_dem_path: str,
-    shpfile: str,
     out_dem_path: str,
-    coreg_method: Coreg | None = NuthKaab() + BiasCorr(bias_func=np.nanmedian),
+    shpfile: str | None,
+    coreg_method: Coreg | None = None,
+    hmode: str = "nuth_kaab",
+    vmode: str = "median",
+    deramp_degree: int = 1,
     grid: str = "ref",
     filtering: bool = True,
     slope_lim: list(AnyNumber) = (0.1, 40),
@@ -2002,9 +2017,12 @@ def dem_coregistration(
 
     :param src_dem_path: path to the input DEM to be coregistered
     :param ref_dem: path to the reference DEM
-    :param shpfile: path to a vector file containing areas to be masked for coregistration
     :param out_dem_path: Path where to save the coregistered DEM
+    :param shpfile: path to a vector file containing areas to be masked for coregistration
     :param coreg_method: The xdem coregistration method, or pipeline. If set to None, DEMs will be resampled to ref grid and optionally filtered, but not coregistered.
+    :param hmode: The method to be used for horizontally aligning the DEMs, e.g. Nuth & Kaab or ICP. Can be any of {list(vmodes_dict.keys())}.
+    :param vmode: The method to be used for vertically aligning the DEMs, e.g. mean/median bias correction or deramping. Can be any of {list(hmodes_dict.keys())}.
+    :param deramp_degree: The degree of the polynomial for deramping.
     :param grid: the grid to be used during coregistration, set either to "ref" or "src".
     :param filtering: if set to True, filtering will be applied prior to coregistration
     :param plot: Set to True to plot a figure of elevation diff before/after coregistration
@@ -2013,7 +2031,19 @@ def dem_coregistration(
 
     :returns: a tuple containing - basename of coregistered DEM, [count of obs, median and NMAD over stable terrain, coverage over roi] before coreg, [same stats] after coreg
     """
+    # Check input arguments
+    if (coreg_method is not None) and ((hmode is not None) or (vmode is not None)):
+        warnings.warn("Both `coreg_method` and `hmode/vmode` are set. Using coreg_method.")
+
+    if hmode not in list(hmodes_dict.keys()):
+        raise ValueError(f"vhmode must be in {list(hmodes_dict.keys())}")
+
+    if vmode not in list(vmodes_dict.keys()):
+        raise ValueError(f"vmode must be in {list(vmodes_dict.keys())}")
+
     # Load both DEMs
+    if verbose:
+        print("Loading and reprojecting input data")
     ref_dem = xdem.DEM(ref_dem_path)
     src_dem = xdem.DEM(src_dem_path)
 
@@ -2026,8 +2056,11 @@ def dem_coregistration(
         raise ValueError(f"`grid` must be either 'ref' or 'src' - currently set to {grid}")
 
     # Create raster mask
-    outlines = gu.Vector(shpfile)
-    stable_mask = ~outlines.create_mask(src_dem)
+    if shpfile is not None:
+        outlines = gu.Vector(shpfile)
+        stable_mask = ~outlines.create_mask(src_dem)
+    else:
+        stable_mask = np.ones(src_dem.data.shape, dtype="bool")
 
     # Calculate dDEM
     ddem = src_dem - ref_dem
@@ -2047,7 +2080,8 @@ def dem_coregistration(
 
     # Calculate dDEM statistics on pixels used for coreg
     inlier_data = ddem.data[inlier_mask].compressed()
-    nstable_orig, med_orig, nmad_orig = len(inlier_data), np.median(inlier_data), xdem.spatialstats.nmad(inlier_data)
+    nstable_orig, mean_orig = len(inlier_data), np.mean(inlier_data)
+    med_orig, nmad_orig = np.median(inlier_data), xdem.spatialstats.nmad(inlier_data)
 
     # Coregister to reference - Note: this will spread NaN
     # Better strategy: calculate shift, update transform, resample
@@ -2055,12 +2089,24 @@ def dem_coregistration(
         coreg_method.fit(ref_dem, src_dem, inlier_mask, verbose=verbose)
         dem_coreg = coreg_method.apply(src_dem, dilate_mask=False)
     elif coreg_method is None:
-        dem_coreg = src_dem
+        # Horizontal coregistration
+        hcoreg_method = hmodes_dict[hmode]
+        hcoreg_method.fit(ref_dem, src_dem, inlier_mask, verbose=verbose)
+        dem_hcoreg = hcoreg_method.apply(src_dem, dilate_mask=False)
+
+        # Vertical coregistration
+        vcoreg_method = vmodes_dict[vmode]
+        if vmode == 'deramp':
+            vcoreg_method.degree = deramp_degree
+        vcoreg_method.fit(ref_dem, dem_hcoreg, inlier_mask, verbose=verbose)
+        dem_coreg = vcoreg_method.apply(dem_hcoreg, dilate_mask=False)
+
     ddem_coreg = dem_coreg - ref_dem
 
     # Calculate new stats
     inlier_data = ddem_coreg.data[inlier_mask].compressed()
-    nstable_coreg, med_coreg, nmad_coreg = len(inlier_data), np.median(inlier_data), xdem.spatialstats.nmad(inlier_data)
+    nstable_coreg, mean_coreg = len(inlier_data), np.mean(inlier_data)
+    med_coreg, nmad_coreg = np.median(inlier_data), xdem.spatialstats.nmad(inlier_data)
 
     # Plot results
     if plot:
@@ -2073,13 +2119,13 @@ def dem_coregistration(
         plt.imshow(ddem.data.squeeze(), cmap="coolwarm_r", vmin=-vmax, vmax=vmax)
         cb = plt.colorbar()
         cb.set_label("Elevation change (m)")
-        ax1.set_title(f"Before coreg\n\nmed = {med_orig:.2f} m - NMAD = {nmad_orig:.2f} m")
+        ax1.set_title(f"Before coreg\n\nmean = {mean_orig:.1f} m - med = {med_orig:.1f} m - NMAD = {nmad_orig:.1f} m")
 
         ax2 = plt.subplot(122, sharex=ax1, sharey=ax1)
         plt.imshow(ddem_coreg.data.squeeze(), cmap="coolwarm_r", vmin=-vmax, vmax=vmax)
         cb = plt.colorbar()
         cb.set_label("Elevation change (m)")
-        ax2.set_title(f"After coreg\n\nmed = {med_coreg:.2f} m - NMAD = {nmad_coreg:.2f} m")
+        ax2.set_title(f"After coreg\n\n\nmean = {mean_coreg:.1f} m - med = {med_coreg:.1f} m - NMAD = {nmad_coreg:.1f} m")
 
         plt.tight_layout()
         if out_fig is None:

From 39cf3ac71e3fc31458e08c639e25611db9c2febe Mon Sep 17 00:00:00 2001
From: adehecq <3285905+adehecq@users.noreply.github.com>
Date: Wed, 5 Oct 2022 09:20:46 +0200
Subject: [PATCH 04/11] Minor comment edit

---
 xdem/coreg.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/xdem/coreg.py b/xdem/coreg.py
index 92ff9b2d..d385d0a8 100644
--- a/xdem/coreg.py
+++ b/xdem/coreg.py
@@ -2016,7 +2016,7 @@ def dem_coregistration(
     verbose: bool = False,
 ):
     """
-    Coregister a selected DEM to a reference DEM.
+    A one-line function to coregister a selected DEM to a reference DEM.
     Reads both DEMs, reproject DEM onto ref DEM grid, mask content of shpfile, run the coregistration and save the coregistered DEM as well as some optional figures and returns some statistics.
 
     :param src_dem_path: path to the input DEM to be coregistered
@@ -2046,6 +2046,7 @@ def dem_coregistration(
         raise ValueError(f"vmode must be in {list(vmodes_dict.keys())}")
 
     # Load both DEMs
+    # TODO: load data only in area of overlap
     if verbose:
         print("Loading and reprojecting input data")
     ref_dem = xdem.DEM(ref_dem_path)

From 1a74d16643a9f22b6cbd8560c4bd82adc1921bd0 Mon Sep 17 00:00:00 2001
From: adehecq <3285905+adehecq@users.noreply.github.com>
Date: Thu, 6 Oct 2022 17:58:05 +0200
Subject: [PATCH 05/11] Fix issue with edge artifacts brought up in #314

---
 xdem/coreg.py | 37 ++++++++++++++++++++-----------------
 1 file changed, 20 insertions(+), 17 deletions(-)

diff --git a/xdem/coreg.py b/xdem/coreg.py
index aed89476..c4c7eb94 100644
--- a/xdem/coreg.py
+++ b/xdem/coreg.py
@@ -1501,7 +1501,7 @@ def apply_matrix(
     :param invert: Invert the transformation matrix.
     :param centroid: The X/Y/Z transformation centroid. Irrelevant for pure translations. Defaults to the midpoint (Z=0)
     :param resampling: The resampling method to use. Can be `nearest`, `bilinear`, `cubic` or an integer from 0-5.
-    :param dilate_mask: Dilate the nan mask to exclude edge pixels that could be wrong.
+    :param dilate_mask: DEPRECATED - This option does not do anything anymore. Will be removed in the future.
 
     :returns: The transformed DEM with NaNs as nodata values (replaces a potential mask of the input `dem`).
     """
@@ -1577,9 +1577,11 @@ def apply_matrix(
     # Shift the elevation values of the soon-to-be-warped DEM.
     filled_dem -= deramp(x_coords, y_coords)
 
-    # Create gap-free arrays of x and y coordinates to be converted into index coordinates.
-    x_inds = rio.fill.fillnodata(transformed_points[:, :, 0].copy(), mask=(~nan_mask).astype("uint8"))
-    y_inds = rio.fill.fillnodata(transformed_points[:, :, 1].copy(), mask=(~nan_mask).astype("uint8"))
+    # Create arrays of x and y coordinates to be converted into index coordinates.
+    x_inds = transformed_points[:, :, 0].copy()
+    x_inds[x_inds == 0] = np.nan
+    y_inds = transformed_points[:, :, 1].copy()
+    y_inds[y_inds == 0] = np.nan
 
     # Divide the coordinates by the resolution to create index coordinates.
     x_inds /= resolution
@@ -1599,19 +1601,20 @@ def apply_matrix(
         transformed_dem = skimage.transform.warp(
             filled_dem, inds, order=resampling_order, mode="constant", cval=np.nan, preserve_range=True
         )
-        # Warp the NaN mask, setting true to all values outside the new frame.
-        tr_nan_mask = (
-            skimage.transform.warp(
-                nan_mask.astype("uint8"), inds, order=resampling_order, mode="constant", cval=1, preserve_range=True
-            )
-            > 0
-        )
-
-    if dilate_mask:
-        tr_nan_mask = scipy.ndimage.binary_dilation(tr_nan_mask, iterations=resampling_order)
-
-    # Apply the transformed nan_mask
-    transformed_dem[tr_nan_mask] = np.nan
+    # TODO: remove these lines when dilate_mask is deprecated
+    #    # Warp the NaN mask, setting true to all values outside the new frame.
+    #     tr_nan_mask = (
+    #         skimage.transform.warp(
+    #             nan_mask.astype("uint8"), inds, order=resampling_order, mode="constant", cval=1, preserve_range=True
+    #         )
+    #         > 0
+    #     )
+
+    # if dilate_mask:
+    #     tr_nan_mask = scipy.ndimage.binary_dilation(tr_nan_mask, iterations=resampling_order)
+
+    # # Apply the transformed nan_mask
+    # transformed_dem[tr_nan_mask] = np.nan
 
     assert np.count_nonzero(~np.isnan(transformed_dem)) > 0, "Transformed DEM has all nans."
 

From 22d3425c1349827e3389f61acddab51395b0fa6d Mon Sep 17 00:00:00 2001
From: adehecq <3285905+adehecq@users.noreply.github.com>
Date: Thu, 6 Oct 2022 18:12:12 +0200
Subject: [PATCH 06/11] Add an experimatal function to slightly extrapolate DEM
 prior to applying transformation

---
 xdem/coreg.py | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/xdem/coreg.py b/xdem/coreg.py
index c4c7eb94..96f7f8a7 100644
--- a/xdem/coreg.py
+++ b/xdem/coreg.py
@@ -1481,6 +1481,7 @@ def apply_matrix(
     centroid: tuple[float, float, float] | None = None,
     resampling: int | str = "bilinear",
     dilate_mask: bool = False,
+    fill_max_search: int = 0,
 ) -> NDArrayf:
     """
     Apply a 3D transformation matrix to a 2.5D DEM.
@@ -1502,6 +1503,9 @@ def apply_matrix(
     :param centroid: The X/Y/Z transformation centroid. Irrelevant for pure translations. Defaults to the midpoint (Z=0)
     :param resampling: The resampling method to use. Can be `nearest`, `bilinear`, `cubic` or an integer from 0-5.
     :param dilate_mask: DEPRECATED - This option does not do anything anymore. Will be removed in the future.
+    :param fill_max_search: Set to > 0 value to fill the DEM before applying the transformation, to avoid spreading\
+    gaps. The DEM will be filled with rasterio.fill.fillnodata with max_search_distance set to fill_max_search.\
+    This is experimental, use at your own risk !
 
     :returns: The transformed DEM with NaNs as nodata values (replaces a potential mask of the input `dem`).
     """
@@ -1534,8 +1538,11 @@ def apply_matrix(
 
     nan_mask = spatial_tools.get_mask(dem)
     assert np.count_nonzero(~nan_mask) > 0, "Given DEM had all nans."
-    # Create a filled version of the DEM. (skimage doesn't like nans)
-    filled_dem = np.where(~nan_mask, demc, np.nan)
+    # Optionally, fill DEM around gaps to reduce spread of gaps
+    if fill_max_search > 0:
+        filled_dem = rio.fill.fillnodata(demc, mask=(~nan_mask).astype("uint8"), max_search_distance=fill_max_search)
+    else:
+        filled_dem = demc  # np.where(~nan_mask, demc, np.nan)  # I don't know why this was needed - to delete
 
     # Get the centre coordinates of the DEM pixels.
     x_coords, y_coords = _get_x_and_y_coords(demc.shape, transform)

From 961cdfc3af8587a3026006e46b13f7f3c865bcb8 Mon Sep 17 00:00:00 2001
From: adehecq <3285905+adehecq@users.noreply.github.com>
Date: Thu, 13 Oct 2022 14:26:57 +0200
Subject: [PATCH 07/11] Load DEMs only in area of overlap using new geoutils
 functionalities.

---
 xdem/coreg.py | 19 +++++++++----------
 1 file changed, 9 insertions(+), 10 deletions(-)

diff --git a/xdem/coreg.py b/xdem/coreg.py
index 96f7f8a7..e84ca456 100644
--- a/xdem/coreg.py
+++ b/xdem/coreg.py
@@ -2167,8 +2167,8 @@ def warp_dem(
 def dem_coregistration(
     src_dem_path: str,
     ref_dem_path: str,
-    out_dem_path: str,
-    shpfile: str | None,
+    out_dem_path: str | None = None,
+    shpfile: str | None = None,
     coreg_method: Coreg | None = None,
     hmode: str = "nuth_kaab",
     vmode: str = "median",
@@ -2186,7 +2186,7 @@ def dem_coregistration(
 
     :param src_dem_path: path to the input DEM to be coregistered
     :param ref_dem: path to the reference DEM
-    :param out_dem_path: Path where to save the coregistered DEM
+    :param out_dem_path: Path where to save the coregistered DEM. If set to None (default), will not save to file.
     :param shpfile: path to a vector file containing areas to be masked for coregistration
     :param coreg_method: The xdem coregistration method, or pipeline. If set to None, DEMs will be resampled to ref grid and optionally filtered, but not coregistered.
     :param hmode: The method to be used for horizontally aligning the DEMs, e.g. Nuth & Kaab or ICP. Can be any of {list(vmodes_dict.keys())}.
@@ -2211,20 +2211,19 @@ def dem_coregistration(
         raise ValueError(f"vmode must be in {list(vmodes_dict.keys())}")
 
     # Load both DEMs
-    # TODO: load data only in area of overlap
     if verbose:
         print("Loading and reprojecting input data")
-    ref_dem = xdem.DEM(ref_dem_path)
-    src_dem = xdem.DEM(src_dem_path)
-
-    # Reproject to common grid
     if grid == "ref":
-        src_dem = src_dem.reproject(ref_dem, resampling="bilinear", silent=True)
+        ref_dem, src_dem = gu.spatial_tools.load_multiple_rasters([ref_dem_path, src_dem_path], ref_grid=0)
     elif grid == "src":
-        ref_dem = ref_dem.reproject(src_dem, resampling="bilinear", silent=True)
+        ref_dem, src_dem = gu.spatial_tools.load_multiple_rasters([ref_dem_path, src_dem_path], ref_grid=1)
     else:
         raise ValueError(f"`grid` must be either 'ref' or 'src' - currently set to {grid}")
 
+    # Convert to DEM instance with Float32 dtype
+    ref_dem = xdem.DEM(ref_dem.astype(np.float32))
+    src_dem = xdem.DEM(src_dem.astype(np.float32))
+
     # Create raster mask
     if shpfile is not None:
         outlines = gu.Vector(shpfile)

From 41b1d91e6ca5321efe8cf9afbcb2ff216daaf76e Mon Sep 17 00:00:00 2001
From: adehecq <3285905+adehecq@users.noreply.github.com>
Date: Thu, 13 Oct 2022 14:41:21 +0200
Subject: [PATCH 08/11] Linting and mypy

---
 xdem/coreg.py | 19 ++++++++++++-------
 1 file changed, 12 insertions(+), 7 deletions(-)

diff --git a/xdem/coreg.py b/xdem/coreg.py
index e84ca456..7ff6ef46 100644
--- a/xdem/coreg.py
+++ b/xdem/coreg.py
@@ -2175,22 +2175,26 @@ def dem_coregistration(
     deramp_degree: int = 1,
     grid: str = "ref",
     filtering: bool = True,
-    slope_lim: list(AnyNumber) = (0.1, 40),
+    slope_lim: list[AnyNumber] | tuple[AnyNumber, AnyNumber] = (0.1, 40),
     plot: bool = False,
     out_fig: str = None,
     verbose: bool = False,
-):
+) -> tuple[xdem.DEM, pd.DataFrame]:
     """
     A one-line function to coregister a selected DEM to a reference DEM.
-    Reads both DEMs, reproject DEM onto ref DEM grid, mask content of shpfile, run the coregistration and save the coregistered DEM as well as some optional figures and returns some statistics.
+    Reads both DEMs, reproject DEM onto ref DEM grid, mask content of shpfile, run the coregistration and save \
+the coregistered DEM as well as some optional figures and returns some statistics.
 
     :param src_dem_path: path to the input DEM to be coregistered
     :param ref_dem: path to the reference DEM
     :param out_dem_path: Path where to save the coregistered DEM. If set to None (default), will not save to file.
     :param shpfile: path to a vector file containing areas to be masked for coregistration
-    :param coreg_method: The xdem coregistration method, or pipeline. If set to None, DEMs will be resampled to ref grid and optionally filtered, but not coregistered.
-    :param hmode: The method to be used for horizontally aligning the DEMs, e.g. Nuth & Kaab or ICP. Can be any of {list(vmodes_dict.keys())}.
-    :param vmode: The method to be used for vertically aligning the DEMs, e.g. mean/median bias correction or deramping. Can be any of {list(hmodes_dict.keys())}.
+    :param coreg_method: The xdem coregistration method, or pipeline. If set to None, DEMs will be resampled to \
+ref grid and optionally filtered, but not coregistered.
+    :param hmode: The method to be used for horizontally aligning the DEMs, e.g. Nuth & Kaab or ICP. Can be any \
+of {list(vmodes_dict.keys())}.
+    :param vmode: The method to be used for vertically aligning the DEMs, e.g. mean/median bias correction or \
+deramping. Can be any of {list(hmodes_dict.keys())}.
     :param deramp_degree: The degree of the polynomial for deramping.
     :param grid: the grid to be used during coregistration, set either to "ref" or "src".
     :param filtering: if set to True, filtering will be applied prior to coregistration
@@ -2198,7 +2202,8 @@ def dem_coregistration(
     :param out_fig: Path to the output figure. If None will display to screen.
     :param verbose: set to True to print details on screen during coregistration.
 
-    :returns: a tuple containing - basename of coregistered DEM, [count of obs, median and NMAD over stable terrain, coverage over roi] before coreg, [same stats] after coreg
+    :returns: a tuple containing - basename of coregistered DEM, [count of obs, median and NMAD over stable \
+terrain, coverage over roi] before coreg, [same stats] after coreg
     """
     # Check input arguments
     if (coreg_method is not None) and ((hmode is not None) or (vmode is not None)):

From c66470b25840f8c68844dd63783976f52b3e1efa Mon Sep 17 00:00:00 2001
From: adehecq <3285905+adehecq@users.noreply.github.com>
Date: Thu, 13 Oct 2022 14:43:29 +0200
Subject: [PATCH 09/11] Upgrade pre-commit

---
 .pre-commit-config.yaml | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 1a7357bc..5b2ab1c8 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -29,7 +29,7 @@ repos:
 
         # Format the code aggressively using black
         - repo: https://github.com/psf/black
-          rev: 22.8.0
+          rev: 22.10.0
           hooks:
                   - id: black
                     args: [--line-length=120]
@@ -52,7 +52,7 @@ repos:
 
         # Lint the code using mypy
         - repo: https://github.com/pre-commit/mirrors-mypy
-          rev: v0.971
+          rev: v0.982
           hooks:
                 - id: mypy
                   args: [
@@ -82,7 +82,7 @@ repos:
 
         # Automatically upgrade syntax to a minimum version
         - repo: https://github.com/asottile/pyupgrade
-          rev: v2.38.0
+          rev: v3.1.0
           hooks:
                 - id: pyupgrade
                   args: [--py37-plus]
@@ -108,6 +108,6 @@ repos:
 
         # Add custom regex lints (see .relint.yml)
         - repo: https://github.com/codingjoe/relint
-          rev: 1.4.0
+          rev: 2.0.0
           hooks:
                 - id: relint

From 79f2e8ac2b50d02fda6243211dfda8a47d387677 Mon Sep 17 00:00:00 2001
From: adehecq <3285905+adehecq@users.noreply.github.com>
Date: Thu, 13 Oct 2022 15:20:48 +0200
Subject: [PATCH 10/11] Improve docstring. Make output file optional.

---
 xdem/coreg.py | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/xdem/coreg.py b/xdem/coreg.py
index 7ff6ef46..5ead1e84 100644
--- a/xdem/coreg.py
+++ b/xdem/coreg.py
@@ -2182,15 +2182,16 @@ def dem_coregistration(
 ) -> tuple[xdem.DEM, pd.DataFrame]:
     """
     A one-line function to coregister a selected DEM to a reference DEM.
-    Reads both DEMs, reproject DEM onto ref DEM grid, mask content of shpfile, run the coregistration and save \
-the coregistered DEM as well as some optional figures and returns some statistics.
+    Reads both DEMs, reprojects them on the same grid, mask content of shpfile, filter steep slopes and outliers, \
+run the coregistration, returns the coregistered DEM and some statistics.
+    Optionally, save the coregistered DEM to file and make a figure.
 
     :param src_dem_path: path to the input DEM to be coregistered
     :param ref_dem: path to the reference DEM
     :param out_dem_path: Path where to save the coregistered DEM. If set to None (default), will not save to file.
     :param shpfile: path to a vector file containing areas to be masked for coregistration
     :param coreg_method: The xdem coregistration method, or pipeline. If set to None, DEMs will be resampled to \
-ref grid and optionally filtered, but not coregistered.
+ref grid and optionally filtered, but not coregistered. Will be used in priority over hmode and vmode.
     :param hmode: The method to be used for horizontally aligning the DEMs, e.g. Nuth & Kaab or ICP. Can be any \
 of {list(vmodes_dict.keys())}.
     :param vmode: The method to be used for vertically aligning the DEMs, e.g. mean/median bias correction or \
@@ -2202,8 +2203,8 @@ def dem_coregistration(
     :param out_fig: Path to the output figure. If None will display to screen.
     :param verbose: set to True to print details on screen during coregistration.
 
-    :returns: a tuple containing - basename of coregistered DEM, [count of obs, median and NMAD over stable \
-terrain, coverage over roi] before coreg, [same stats] after coreg
+    :returns: a tuple containing 1) coregistered DEM as an xdem.DEM instance and 2) DataFrame of coregistration \
+statistics (count of obs, median and NMAD over stable terrain) before and after coreg.
     """
     # Check input arguments
     if (coreg_method is not None) and ((hmode is not None) or (vmode is not None)):
@@ -2311,7 +2312,8 @@ def dem_coregistration(
             plt.close()
 
     # Save coregistered DEM
-    dem_coreg.save(out_dem_path, tiled=True)
+    if out_dem_path is not None:
+        dem_coreg.save(out_dem_path, tiled=True)
 
     # Save stats to DataFrame
     out_stats = pd.DataFrame(

From f4b8177ecf76564d4dd62695fe190047c1cf84e1 Mon Sep 17 00:00:00 2001
From: adehecq <3285905+adehecq@users.noreply.github.com>
Date: Fri, 14 Oct 2022 16:28:37 +0200
Subject: [PATCH 11/11] Add a 5 NMAD filter to the filtering option.

---
 xdem/coreg.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/xdem/coreg.py b/xdem/coreg.py
index 5ead1e84..94322416 100644
--- a/xdem/coreg.py
+++ b/xdem/coreg.py
@@ -2242,8 +2242,10 @@ def dem_coregistration(
 
     # Filter gross outliers in stable terrain
     if filtering:
-        # TO DO implement the NMAD filter in xdem
-        inlier_mask = stable_mask  # nmad_filter(ddem.data, stable_mask, verbose=False)
+        # Remove gross blunders where dh differ by 5 NMAD from the median
+        inlier_mask = stable_mask & (np.abs(ddem.data - np.median(ddem)) < 5 * xdem.spatialstats.nmad(ddem)).filled(
+            False
+        )
 
         # Exclude steep slopes for coreg
         slope = xdem.terrain.slope(ref_dem)