Fixed issue #56 (Setting align_grids to False in case of global...

Fixed issue #56

 (Setting align_grids to False in case of global co-registration does not avoid resampling.).

Signed-off-by: Daniel Scheffler <danschef@gfz-potsdam.de>

arosics/CoReg.py

@@ -184,10 +184,12 @@ class COREG(object):
         """Get an instance of the COREG class.
 
         :param im_ref:
-            source path (any GDAL compatible image format is supported) or GeoArray instance of reference image
+            source path or GeoArray instance of reference image
+            (any GDAL compatible image format is supported)
 
         :param im_tgt:
-            source path (any GDAL compatible image format is supported) or GeoArray instance of image to be shifted
+            source path or GeoArray instance of the target image, i.e., the image to be shifted
+            (any GDAL compatible image format is supported)
 
         :param path_out:
             target path of the coregistered image
@@ -222,14 +224,21 @@ class COREG(object):
             maximum shift distance in reference image pixel units (default: 5 px)
 
         :param align_grids:
-            align the coordinate grids of the image to be and the reference image (default: False)
+            True: align the input coordinate grid to the reference (does not affect the output pixel size as long as
+            input and output pixel sizes are compatible (5:30 or 10:30 but not 4:30), default = False
+
+            - NOTE: If this is set to False, the mis-registration in the target image is corrected by updating its
+              geocoding information only, i.e., without performing any resampling which preserves the original
+              image values (except that a resampling is needed due to other reasons, e.g., if 'match_gsd',
+              'out_gsd' or 'target_xyGrid' are given).
 
         :param match_gsd:
-            match the output pixel size to pixel size of the reference image (default: False)
+            True: match the input pixel size to the reference pixel size,
+            default = False
 
         :param out_gsd:
-            xgsd ygsd: set the output pixel size in map units
-            (default: original pixel size of the image to be shifted)
+            output pixel size (X/Y) in units of the reference coordinate system (default = pixel size of the target
+            image), given values are overridden by match_gsd=True
 
         :param target_xyGrid:
             a list with a target x-grid and a target y-grid like [[15,45], [15,45]]
@@ -247,7 +256,7 @@ class COREG(object):
 
         :param footprint_poly_ref:
             footprint polygon of the reference image (WKT string or shapely.geometry.Polygon),
-             e.g. 'POLYGON ((299999 6000000, 299999 5890200, 409799 5890200, 409799 6000000, 299999 6000000))'
+            e.g. 'POLYGON ((299999 6000000, 299999 5890200, 409799 5890200, 409799 6000000, 299999 6000000))'
 
         :param footprint_poly_tgt:
             footprint polygon of the image to be shifted (WKT string or shapely.geometry.Polygon)

arosics/DeShifter.py

@@ -173,11 +173,11 @@ class DESHIFTER(object):
             coreg_results['updated map info'] = coreg_results['updated map info means']
 
         # in case of global shift correction -> the updated map info from coreg_results already has the final map info
-        # BUT: this will updated in correct_shifts() if clipextent is given or warping is needed
+        # BUT: this will be updated in correct_shifts() if clipextent is given or warping is needed
         if not self.GCPList:
             mapI = coreg_results['updated map info']
-            self.updated_map_info = mapI if mapI else geotransform2mapinfo(self.shift_gt, self.shift_prj)
-            self.updated_gt = mapinfo2geotransform(self.updated_map_info) if mapI else self.shift_gt
+            self.updated_map_info = mapI or geotransform2mapinfo(self.shift_gt, self.shift_prj)
+            self.updated_gt = mapinfo2geotransform(self.updated_map_info) or self.shift_gt
             self.original_map_info = coreg_results['original map info']
         self.updated_projection = self.ref_prj
 
@@ -249,8 +249,14 @@ class DESHIFTER(object):
                 # use origin of reference image and gsd of input image
                 out_grid = get_grid(self.ref_gt, self.im2shift.xgsd, self.im2shift.ygsd)
             else:
-                # use input image grid
-                out_grid = get_grid(self.im2shift.geotransform, self.im2shift.xgsd, self.im2shift.ygsd)
+                if not self.GCPList:
+                    # in case of global co-registration:
+                    # -> use the target image grid but update the origin (shift-correction without resampling)
+                    out_grid = get_grid(self.updated_gt, self.im2shift.xgsd, self.im2shift.ygsd)
+                else:
+                    # in case of local co-registration:
+                    # -> use input image grid
+                    out_grid = get_grid(self.im2shift.geotransform, self.im2shift.xgsd, self.im2shift.ygsd)
 
         return out_grid
 
@@ -323,6 +329,11 @@ class DESHIFTER(object):
             self.is_resampled = False
             xmin, ymin, xmax, ymax = self._get_out_extent()
 
+            if not self.q:
+                print("NOTE: The detected shift is corrected by updating the map info of the target image only, i.e., "
+                      "without any resampling. Set the 'align_grids' parameter to True if you need the target and the "
+                      "reference coordinate grids to be aligned.")
+
             if self.cliptoextent:
                 # TODO validate results
                 # TODO -> output extent does not seem to be the requested one! (only relevant if align_grids=False)

tests/test_COREG.py

@@ -307,6 +307,30 @@ class CompleteWorkflow_INTER1_S2A_S2A(unittest.TestCase):
             CR.show_matchWin(interactive=True, after_correction=False)
             CR.show_image_footprints()
 
+    def test_correct_shifts_without_resampling(self):
+        kw = self.coreg_kwargs.copy()
+        kw['align_grids'] = False  # =default
+        kw['progress'] = True
+
+        CR = self.run_shift_detection_correction(self.ref_path, self.tgt_path, **kw)
+        self.assertTrue(CR.success)
+        self.assertTrue(CR.deshift_results['is shifted'])
+        self.assertFalse(CR.deshift_results['is resampled'])
+        self.assertTrue(np.array_equal(CR.shift[:], CR.deshift_results['arr_shifted']))
+        self.assertFalse(np.array_equal(np.array(CR.shift.gt), np.array(CR.deshift_results['updated geotransform'])))
+
+    def test_correct_shifts_with_resampling(self):
+        kw = self.coreg_kwargs.copy()
+        kw['align_grids'] = True
+        kw['progress'] = True
+
+        CR = self.run_shift_detection_correction(self.ref_path, self.tgt_path, **kw)
+        self.assertTrue(CR.success)
+        self.assertTrue(CR.deshift_results['is shifted'])
+        self.assertTrue(CR.deshift_results['is resampled'])
+        self.assertFalse(np.array_equal(CR.shift[:], CR.deshift_results['arr_shifted']))
+        self.assertTrue(np.array_equal(np.array(CR.shift.gt), np.array(CR.deshift_results['updated geotransform'])))
+
 
 if __name__ == '__main__':
     import nose2