From 609855a1c47b883d0457f127ff3a868c3d080641 Mon Sep 17 00:00:00 2001 From: Romain Hugonnet Date: Thu, 28 Mar 2024 14:14:15 -0800 Subject: [PATCH] Amaurys comments --- xdem/coreg/affine.py | 1 - xdem/coreg/base.py | 35 ++++++++++++++++++++++------------- xdem/coreg/biascorr.py | 7 ++++--- 3 files changed, 26 insertions(+), 17 deletions(-) diff --git a/xdem/coreg/affine.py b/xdem/coreg/affine.py index a9b29c9b..49747616 100644 --- a/xdem/coreg/affine.py +++ b/xdem/coreg/affine.py @@ -583,7 +583,6 @@ def _fit_rst_pts( ) # If raster was reference, invert the matrix - # TODO: Move matrix/invert_matrix to affine module? if ref == "raster": matrix = xdem.coreg.base.invert_matrix(matrix) diff --git a/xdem/coreg/base.py b/xdem/coreg/base.py index ee532afa..39c34972 100644 --- a/xdem/coreg/base.py +++ b/xdem/coreg/base.py @@ -429,10 +429,10 @@ def _preprocess_coreg_fit_point_point( ) -> tuple[gpd.GeoDataFrame, gpd.GeoDataFrame]: """Pre-processing and checks of fit for point-point input.""" - ref_dem = reference_elev - tba_dem = to_be_aligned_elev.to_crs(crs=reference_elev.crs) + ref_elev = reference_elev + tba_elev = to_be_aligned_elev.to_crs(crs=reference_elev.crs) - return ref_dem, tba_dem + return ref_elev, tba_elev def _preprocess_coreg_fit( @@ -472,16 +472,16 @@ def _preprocess_coreg_fit( point_elev = reference_elev ref = "point" - rst_elev, point_elev, inlier_mask, transform, crs = _preprocess_coreg_fit_raster_point( + raster_elev, point_elev, inlier_mask, transform, crs = _preprocess_coreg_fit_raster_point( raster_elev=raster_elev, point_elev=point_elev, inlier_mask=inlier_mask, transform=transform, crs=crs ) if ref == "raster": - ref_elev = rst_elev + ref_elev = raster_elev tba_elev = point_elev else: ref_elev = point_elev - tba_elev = rst_elev + tba_elev = raster_elev # If both inputs are points, simply reproject to the same CRS else: @@ -555,7 +555,12 @@ def _postprocess_coreg_apply_rst( resample: bool, resampling: rio.warp.Resampling | None = None, ) -> tuple[NDArrayf | gu.Raster, affine.Affine]: - """Post-processing and checks of apply for raster input.""" + """ + Post-processing and checks of apply for raster input. + + Here, "elev" and "transform" corresponds to user input, and are required to transform back the output that is + composed of "applied_elev" and "out_transform". + """ # Ensure the dtype is OK applied_elev = applied_elev.astype("float32") @@ -605,7 +610,12 @@ def _postprocess_coreg_apply( resample: bool, resampling: rio.warp.Resampling | None = None, ) -> tuple[NDArrayf | gpd.GeoDataFrame, affine.Affine]: - """Post-processing and checks of apply for any input.""" + """ + Post-processing and checks of apply for any input. + + Here, "elev" and "transform" corresponds to user input, and are required to transform back the output that is + composed of "applied_elev" and "out_transform". + """ # Define resampling resampling = resampling if isinstance(resampling, rio.warp.Resampling) else resampling_method_from_str(resampling) @@ -1181,9 +1191,6 @@ def fit( if self._meta["subsample"] != 1: self._meta["random_state"] = random_state - # TODO: Add preproc for points too - # TODO: Rename into "checks", because not much is preprocessed in the end - # (has to happen in the _fit_func itself, whether for subsampling or # Pre-process the inputs, by reprojecting and converting to arrays ref_elev, tba_elev, inlier_mask, transform, crs = _preprocess_coreg_fit( reference_elev=reference_elev, @@ -2037,6 +2044,9 @@ def fit( **kwargs: Any, ) -> CoregType: + if isinstance(reference_elev, gpd.GeoDataFrame) and isinstance(to_be_aligned_elev, gpd.GeoDataFrame): + raise NotImplementedError("Blockwise coregistration does not yet support two elevation point cloud inputs.") + # Check if subsample arguments are different from their default value for any of the coreg steps: # get default value in argument spec and "subsample" stored in meta, and compare both are consistent if not isinstance(self.procstep, CoregPipeline): @@ -2065,8 +2075,7 @@ def fit( crs=crs, ) - # TODO: Blockwise can only work if one of the two is a Raster... or by defining a grid somehow? - groups = self.subdivide_array(tba_dem.shape) + groups = self.subdivide_array(tba_dem.shape if isinstance(tba_dem, np.ndarray) else ref_dem.shape) indices = np.unique(groups) diff --git a/xdem/coreg/biascorr.py b/xdem/coreg/biascorr.py index d6cc2a6a..b9e87674 100644 --- a/xdem/coreg/biascorr.py +++ b/xdem/coreg/biascorr.py @@ -392,6 +392,7 @@ def _fit_rst_pts( # type: ignore inlier_pts_alltrue = np.ones(len(pts), dtype=bool) # Below, we derive 1D arrays for the rst_rst function to take over after interpolating to the point coordinates + # (as rst_rst works for 1D arrays as well as 2D arrays, as long as coordinates match) # Convert ref or tba depending on which is the point dataset if isinstance(ref_elev, gpd.GeoDataFrame): @@ -413,7 +414,7 @@ def _fit_rst_pts( # type: ignore else: bias_vars_pts = None - # Send to raster-raster fit + # Send to raster-raster fit but using 1D arrays instead of 2D arrays (flattened anyway during analysis) self._fit_biascorr( ref_elev=ref_elev_pts, tba_elev=tba_elev_pts, @@ -1073,7 +1074,7 @@ def _fit_rst_rst( # type: ignore ) -> None: # The number of parameters in the first guess defines the polynomial order when calling np.polyval2d - p0 = np.ones(shape=((self._meta["poly_order"] + 1) * (self._meta["poly_order"] + 1))) + p0 = np.ones(shape=((self._meta["poly_order"] + 1) **2)) # Coordinates (we don't need the actual ones, just array coordinates) xx, yy = np.meshgrid(np.arange(0, ref_elev.shape[1]), np.arange(0, ref_elev.shape[0])) @@ -1110,7 +1111,7 @@ def _fit_rst_pts( # type: ignore rast_elev = ref_elev if not isinstance(ref_elev, gpd.GeoDataFrame) else tba_elev # The number of parameters in the first guess defines the polynomial order when calling np.polyval2d - p0 = np.ones(shape=((self._meta["poly_order"] + 1) * (self._meta["poly_order"] + 1))) + p0 = np.ones(shape=((self._meta["poly_order"] + 1) ** 2)) # Coordinates (we don't need the actual ones, just array coordinates) xx, yy = np.meshgrid(np.arange(0, rast_elev.shape[1]), np.arange(0, rast_elev.shape[0]))