Fixed issue # 17 (Coregistration sometimes fails in case of floating point...

Fixed issue # 17 (Coregistration sometimes fails in case of floating point coordinates of the input images.)

.gitlab-ci.yml

@@ -23,6 +23,7 @@ test_arosics:
     - source activate ci_env
     - export GDAL_DATA=/root/miniconda3/envs/ci_env/share/gdal
     - export PYTHONPATH=$PYTHONPATH:/root  # /root <- directory needed later
+    - pip install py_tools_ds>=0.13.7  # FIXME remove as soon as docker runner has been updated
     - make nosetests
     - make docs
   artifacts:

arosics/CoReg.py

@@ -4,7 +4,7 @@ import os
 import time
 import warnings
 from copy import copy
-from typing import Iterable
+from typing import Iterable, Union, Tuple  # noqa F401
 
 # custom
 try:
@@ -698,7 +698,8 @@ class COREG(object):
         # Check, ob match Fenster größer als anderes Fenster
         if not (matchBox.mapPoly.within(otherBox.mapPoly) or matchBox.mapPoly == otherBox.mapPoly):
             # dann für anderes Fenster kleinstes Fenster finden, das match-Fenster umgibt
-            otherBox.boxImYX = get_smallest_boxImYX_that_contains_boxMapYX(matchBox.boxMapYX, otherBox.gt)
+            otherBox.boxImYX = get_smallest_boxImYX_that_contains_boxMapYX(
+                matchBox.boxMapYX, otherBox.gt, tolerance_ndigits=5)  # avoids float coordinate rounding issues
 
         # evtl. kann es sein, dass bei Shift-Fenster-Vergrößerung das shift-Fenster zu groß für den overlap wird
         t_start = time.time()
@@ -707,7 +708,8 @@ class COREG(object):
             xLarger, yLarger = otherBox.is_larger_DimXY(overlapWin.boundsIm)
             matchBox.buffer_imXY(-1 if xLarger else 0, -1 if yLarger else 0)
             previous_area = otherBox.mapPoly.area
-            otherBox.boxImYX = get_smallest_boxImYX_that_contains_boxMapYX(matchBox.boxMapYX, otherBox.gt)
+            otherBox.boxImYX = get_smallest_boxImYX_that_contains_boxMapYX(
+                matchBox.boxMapYX, otherBox.gt, tolerance_ndigits=5)  # avoids float coordinate rounding issues)
 
             if previous_area == otherBox.mapPoly.area or time.time() - t_start > 1.5:
                 # happens e.g in case of a triangular footprint
@@ -825,7 +827,7 @@ class COREG(object):
 
     @staticmethod
     def _shrink_winsize_to_binarySize(win_shape_YX, target_size=None):
-        # type: (tuple, tuple, int , int) -> any
+        # type: (tuple, tuple) -> Union[Tuple[int, int], None]
         """Shrinks a given window size to the closest binary window size (a power of 2) -
         separately for X- and Y-dimension.
 
@@ -968,7 +970,7 @@ class COREG(object):
     def _get_grossly_deshifted_images(self, im0, im1, x_intshift, y_intshift):
         # TODO this is also implemented in GeoArray # this should update ref.win.data and shift.win.data
         # FIXME avoid that matching window gets smaller although shifting it  with the previous win_size would not move
-        # FIXME it into nodata-area
+        #       it into nodata-area
         # get_grossly_deshifted_im0
         old_center_YX = np.array(im0.shape) / 2
         new_center_YX = [old_center_YX[0] + y_intshift, old_center_YX[1] + x_intshift]

requirements.txt

 geoarray>=0.8.0
-py_tools_ds>=0.13.4
+py_tools_ds>=0.13.7
 cmocean
 numpy
 gdal

setup.py

@@ -18,7 +18,7 @@ with open("arosics/version.py") as version_file:
     exec(version_file.read(), version)
 
 requirements = ['numpy', 'gdal', 'shapely', 'scikit-image', 'matplotlib', 'geopandas', 'geoarray>=0.8.0',
-                'py_tools_ds>=0.13.4', 'plotly', 'cmocean', 'six', 'folium>=0.6.0', 'geojson'
+                'py_tools_ds>=0.13.7', 'plotly', 'cmocean', 'six', 'folium>=0.6.0', 'geojson'
                 # 'pykrige'  # conda install --yes -c conda-forge pykrige
                 # 'pyfftw', # conda install --yes -c conda-forge pyfftw=0.10.4 ; \
                 # 'basemap', # conda install --yes -c conda-forge basemap; \

tests/CI_docker/context/environment_arosics.yml

@@ -47,4 +47,4 @@ dependencies:
     - cmocean
     - plotly
     - geoarray>=0.8.0
-    - py_tools_ds>=0.13.4
+    - py_tools_ds>=0.13.7

tests/test_COREG.py

@@ -103,6 +103,26 @@ class CompleteWorkflow_INTER1_S2A_S2A(unittest.TestCase):
                                                         footprint_poly_tgt=None))
         self.assertTrue(CR.success)
 
+    def test_shift_calculation_with_float_coords(self):
+        """Test with default parameters - should compute X/Y shifts properly ad write the de-shifted target image."""
+
+        # overwrite gt and prj
+        ref = GeoArray(self.ref_path)
+        ref.to_mem()
+        ref.filePath = None
+        ref.gt = [330000.00000001, 10.1, 0.0, 5862000.0000001, 0.0, -10.1]
+        tgt = GeoArray(self.tgt_path)
+        tgt.to_mem()
+        tgt.filePath = None
+        tgt.gt = [335440.0000001, 10.1, 0.0, 5866490.0000001, 0.0, -10.1]
+
+        CR = self.run_shift_detection_correction(ref, tgt,
+                                                 **dict(self.coreg_kwargs,
+                                                        wp=(341500.0, 5861440.0),
+                                                        footprint_poly_ref=None,
+                                                        footprint_poly_tgt=None))
+        self.assertTrue(CR.success)
+
     def test_shift_calculation_inmem_gAs_path_out_auto(self):
         """Test input parameter path_out='auto' in case input reference/ target image are in-memory GeoArrays."""
         ref = GeoArray(np.random.randint(1, 100, (1000, 1000)))

tests/test_COREG_LOCAL.py

@@ -78,6 +78,28 @@ class CompleteWorkflow_INTER1_S2A_S2A(unittest.TestCase):
         self.assertTrue(os.path.exists(self.coreg_kwargs['path_out']),
                         'Output of local co-registration has not been written.')
 
+    def test_calculation_of_tie_point_grid_float_coords(self):
+        # NOTE: This does not test against unequaly sized output of get_image_windows_to_match().
+
+        # overwrite gt and prj
+        ref = GeoArray(self.ref_path)
+        ref.to_mem()
+        ref.filePath = None
+        ref.gt = [330000.19999996503, 0.6, 0.0, 5862000.7999997628, 0.0, -0.6]
+        # ref.gt = [330000.1, 10.1, 0.0, 5862000.1, 0.0, -10.1]
+        tgt = GeoArray(self.tgt_path)
+        tgt.to_mem()
+        tgt.filePath = None
+        tgt.gt = [330000.19999996503, 0.6, 0.0, 5862000.7999997628, 0.0, -0.6]
+
+        # get instance of COREG_LOCAL object
+        CRL = COREG_LOCAL(ref, tgt, **dict(**self.coreg_kwargs))
+
+        # use the getter of the CoRegPoints_table to calculate tie point grid
+        # noinspection PyStatementEffect
+        CRL.CoRegPoints_table
+
+
 # if __name__ == '__main__':
 #     unittest.main(argv=['first-arg-is-ignored'],exit=False, verbosity=2)
 #