Fix directional bias example

doc/source/biascorr.md

@@ -143,7 +143,7 @@ Several good practices help performing a successful bias correction:
 - **Avoid using "fit" with a subsample size larger than 1,000,000:** Otherwise the optimizer will be extremely slow and might fail with a memory error; consider using "bin_and_fit" instead to reduce the data size before the optimization which still allows to utilize all the data,
 - **Avoid using "fit" or "bin_and_fit" for more than 2 dimensions (input variables):** Fitting a parametric form in more than 2 dimensions is quite delicate, consider using "bin" or sequential 1D corrections instead,
 - **Use at least 1000 bins for all dimensions, being mindful about dimension number:** Using a small bin size is generally too rough, but a large bin size will grow exponentially with the number of bias variables,
-- **Use customized bin edges for data with large extreme values:** Passing simply a bin size will set the min/max of the data as the full binning range, which can be impractical (e.g., most curvatures lies between -2/2 but can take values of 10000+).
+- **Use customized bin edges for data with large extreme values:** Passing simply a bin size will set the min/max of the data as the full binning range, which can be impractical (e.g., most curvatures lie between -2/2 but can take values of 10000+).
 
 ## Bias-correction methods
 

tests/test_coreg/test_biascorr.py

@@ -568,7 +568,7 @@ def test_terrainbias(self) -> None:
         tb = biascorr.TerrainBias()
 
         assert tb.meta["inputs"]["fitorbin"]["fit_or_bin"] == "bin"
-        assert tb.meta["inputs"]["fitorbin"]["bin_sizes"] == 1000
+        assert tb.meta["inputs"]["fitorbin"]["bin_sizes"] == 100
         assert tb.meta["inputs"]["fitorbin"]["bin_statistic"] == np.nanmedian
         assert tb.meta["inputs"]["specific"]["terrain_attribute"] == "maximum_curvature"
         assert tb._needs_vars is False

xdem/coreg/base.py

@@ -91,11 +91,11 @@
     "tolerance": "Tolerance to reach (pixel size)",
     "last_iteration": "Iteration at which algorithm stopped",
     "all_tolerances": "Tolerances at each iteration",
-    "terrain_attribute": "Terrain attribute used for TerrainBias",
-    "angle": "Angle used for DirectionalBias",
-    "poly_order": "Polynomial order used for Deramp",
-    "best_poly_order": "Best polynomial order kept for fit",
-    "best_nb_sin_freq": "Best number of sinusoid frequencies kept for fit",
+    "terrain_attribute": "Terrain attribute used for correction",
+    "angle": "Angle of directional correction",
+    "poly_order": "Polynomial order",
+    "best_poly_order": "Best polynomial order",
+    "best_nb_sin_freq": "Best number of sinusoid frequencies",
     "vshift_reduc_func": "Reduction function used to remove vertical shift",
     "centroid": "Centroid found for affine rotation",
     "shift_x": "Eastward shift estimated (georeferenced unit)",

xdem/coreg/biascorr.py

@@ -12,7 +12,7 @@
 import xdem.spatialstats
 from xdem._typing import NDArrayb, NDArrayf
 from xdem.coreg.base import Coreg, fit_workflows
-from xdem.fit import polynomial_2d
+from xdem.fit import polynomial_2d, sumsin_1d
 
 BiasCorrType = TypeVar("BiasCorrType", bound="BiasCorr")
 
@@ -356,12 +356,6 @@ def _fit_rst_rst(  # type: ignore
             along_track_angle=self._meta["inputs"]["specific"]["angle"],
         )
 
-        # Parameters dependent on resolution cannot be derived from the rotated x coordinates, need to be passed below
-        if "hop_length" not in kwargs:
-            # The hop length will condition jump in function values, need to be larger than average resolution
-            average_res = (transform[0] + abs(transform[4])) / 2
-            kwargs.update({"hop_length": average_res})
-
         super()._fit_rst_rst_and_rst_pts(
             ref_elev=ref_elev,
             tba_elev=tba_elev,

xdem/fit.py

@@ -498,9 +498,10 @@ def wrapper_cost_sumofsin(p: NDArrayf, x: NDArrayf, y: NDArrayf) -> float:
         return _cost_sumofsin(x, y, cost_func, *p)
 
     # If no significant resolution is provided, assume that it is the mean difference between sampled X values
-    if hop_length is None:
-        x_res = np.mean(np.diff(np.sort(xdata)))
-        hop_length = x_res
+    x_res = np.mean(np.diff(np.sort(xdata)))
+
+    # The hop length will condition jump in function values, needs of magnitude slightly lower than the signal
+    hop_length = (np.percentile(ydata, 90) - np.percentile(ydata, 10))
 
     # Loop on all frequencies
     costs = np.empty(max_nb_frequency)
@@ -522,7 +523,7 @@ def wrapper_cost_sumofsin(p: NDArrayf, x: NDArrayf, y: NDArrayf) -> float:
             ub_phase = 2 * np.pi
             # For the wavelength: from the resolution and coordinate extent
             # (we don't want the lower bound to be zero, to avoid divisions by zero)
-            lb_wavelength = hop_length / 5
+            lb_wavelength = x_res / 5
             ub_wavelength = xdata.max() - xdata.min()
 
             b = []
@@ -543,12 +544,12 @@ def wrapper_cost_sumofsin(p: NDArrayf, x: NDArrayf, y: NDArrayf) -> float:
             print(ub)
 
         # Minimize the globalization with a larger number of points
-        minimizer_kwargs = dict(args=(xdata, ydata), method="L-BFGS-B", bounds=scipy_bounds)
+        minimizer_kwargs = dict(args=(xdata, ydata), bounds=scipy_bounds)
         myresults = scipy.optimize.basinhopping(
             wrapper_cost_sumofsin,
             p0,
             disp=verbose,
-            T=hop_length * 50,
+            T=hop_length,
             minimizer_kwargs=minimizer_kwargs,
             seed=random_state,
             **kwargs,