GlacioHack / xdem

Analysis of digital elevation models (DEMs)
https://xdem.readthedocs.io/
Apache License 2.0
148 stars 41 forks source link

[POC]: Add an initial shift for coregistration #603

Open adebardo opened 1 month ago

adebardo commented 1 month ago

Context

The CNES wants to integrate the ability to instantiate an initial offset for the co-registration calculation. Because the dem_coregistration function allows for easy addition of this feature and we will recommend its use to CNES members, we propose simple modifications at this level.

Implementation

  1. Using the translate function from Geoutils:

    The translate function from Geoutils allows for translating a Digital Elevation Model (DEM). We will call this function in the dem_coregsitration method .

def dem_coregistration(
    src_dem_path: str | RasterType,
    ref_dem_path: str | RasterType,
    out_dem_path: str | None = None,
    coreg_method: Coreg | None = NuthKaab() + VerticalShift(),
    grid: str = "ref",
    resample: bool = False,
    resampling: rio.warp.Resampling | None = rio.warp.Resampling.bilinear,
    shp_list: list[str | gu.Vector] | tuple[str | gu.Vector] | tuple[()] = (),
    inout: list[int] | tuple[int] | tuple[()] = (),
    filtering: bool = True,
    dh_max: Number = None,
    nmad_factor: Number = 5,
    slope_lim: list[Number] | tuple[Number, Number] = (0.1, 40),
    plot: bool = False,
    out_fig: str = None,
+   estimated_initial_shift : list = [0, 0] 
) -> tuple[DEM, Coreg, pd.DataFrame, NDArrayf]:
+   """
+   explain estimated_initial_shift + affine only
+   """

        if isinstance(ref_dem_path, str):
            if grid == "ref":
                ref_dem, src_dem = gu.raster.load_multiple_rasters([ref_dem_path, src_dem_path], ref_grid=0)
            elif grid == "src":
                ref_dem, src_dem = gu.raster.load_multiple_rasters([ref_dem_path, src_dem_path], ref_grid=1)
        else:
            ref_dem = ref_dem_path
            src_dem = src_dem_path
            if grid == "ref":
                src_dem = src_dem.reproject(ref_dem, silent=True)
            elif grid == "src":
                ref_dem = ref_dem.reproject(src_dem, silent=True)
+      if estimated_initial_shift != [0, 0] :
+          logging.warning("Initial shift in affine mode only")
+          estimated_initial_shift = [estimated_initial_shift[0] * pixel_resolution_x, estimated_initial_shift[0] * pixel_resolution_y]
+          src_dem.translate(estimated_initial_shift[0], estimated_initial_shift[1], inplace=True)
    .
    .
    .
    # Coregister to reference - Note: this will spread NaN
    coreg_method.fit(ref_dem, src_dem, inlier_mask)
    dem_coreg = coreg_method.apply(src_dem, resample=resample, resampling=resampling)
+   if estimated_initial_shift != [0, 0] :
+       calculated_shift_x = dem_coreg.meta["shift_x"]
+       calculated_shift_y = dem_coreg.meta["shift_y"]

+       dem_coreg.meta["shift_x"] = calculated_shift_x + estimated_initial_shift[0] 
+       dem_coreg.meta["shift_y"] = calculated_shift_x + estimated_initial_shift[1]

Tests

Documentation

/estimation 3d

rhugonnet commented 3 weeks ago

Looks good! :slightly_smiling_face:

So the "explain affine only" would essentially be a check like coreg_method.is_affine for a Coreg object or all(c.is_affine for c in coreg_method) for a CoregPipeline object? I'm not sure of the behaviour of the is_affine property on a pipeline right now (it might fail), we could modify it so it derives the proper attribute directly here: https://github.com/GlacioHack/xdem/blob/595acb9bf2ba150232cf161fc09f8747c84b5602/xdem/coreg/base.py#L1602 (either by finding a solution for it to run on all subclasses, or simply by overridding the property in CoregPipeline).