Merge branch 'master' into feature/add_multi_projection_support

.github/create_release_from_gitlab_ci.sh

@@ -14,6 +14,9 @@ API_JSON=$(printf '{"tag_name":"%s",
 AUTH_HEADER="Authorization: token $GITHUB_RELEASE_TOKEN"
 curl $URL_RELEASES --data "$API_JSON" -H "$AUTH_HEADER"
 
+# sleep 5 seconds, otherwise the check below may fail
+sleep 5
+
 # get latest release tag name
 LATEST_RELEASE=$(curl -s "$URL_RELEASES"/latest | grep -oP '"tag_name": "\K(.*)(?=")')
 

HISTORY.rst

@@ -2,6 +2,26 @@
 History
 =======
 
+1.7.0 (2021-09-30)
+------------------
+
+* Added method Tie_Point_Grid.calc_overall_stats() + test to compute overall statistics from all tie points found.
+
+
+1.6.2 (2021-09-29)
+------------------
+
+* Fixed 'too many values to unpack' exception in COREG_LOCAL.view_CoRegPoints().
+* Added new parameters to Tie_Point_Grid.plot_shift_distribution().
+* Added documentation to COREG_LOCAL.view_CoRegPoints().
+
+
+1.6.1 (2021-09-29)
+------------------
+
+* The output map of COREG_LOCAL.view_CoRegPoints() is now cropped nicely. Added parameter 'figsize_multiplier'.
+
+
 1.6.0 (2021-09-27)
 ------------------
 

arosics/CoReg_local.py

@@ -498,6 +498,7 @@ class COREG_LOCAL(object):
                          backgroundIm: str = 'tgt',
                          hide_filtered: bool = True,
                          figsize: tuple = None,
+                         figsize_multiplier: float = 1,
                          title: str = '',
                          vector_scale: float = 1.,
                          savefigPath: str = '',
@@ -518,36 +519,45 @@ class COREG_LOCAL(object):
         :param backgroundIm:        whether to use the target or the reference image as map background. Possible
                                     options are 'ref' and 'tgt' (default: 'tgt')
         :param hide_filtered:       hide all points that have been filtered out according to tie point filter level
-        :param figsize:             size of the figure to be viewed, e.g. (10, 10)
+        :param figsize:             size of the figure to be viewed, e.g. (10, 10); automatically estimated if not given
+        :param figsize_multiplier:  if given, the figure size is multiplied with this value
         :param title:               plot title
         :param vector_scale:        scale factor for shift vector length (default: 1 -> no scaling)
-        :param savefigPath:
-        :param savefigDPI:
-        :param showFig:             whether to show or to hide the figure
-        :param vmin:
-        :param vmax:
-        :param return_map:
-        :return:
+        :param savefigPath:         path where to save the figure
+        :param savefigDPI:          DPI resolution of the output figure when saved to disk (default: 96)
+        :param showFig:             whether to show or to hide the figure (default: True)
+        :param vmin:                minimum value of 'attribute2plot' to be included in the figure
+        :param vmax:                maximum value of 'attribute2plot' to be included in the figure
+        :param return_map:          whether to return the figure and axis objects (default: False)
+        :return:    tuple of figure and axis objects or None in case return_map is set to False
         """
         from matplotlib import pyplot as plt  # noqa
         from matplotlib.offsetbox import AnchoredText
         from cartopy.crs import PlateCarree
         from mpl_toolkits.axes_grid1 import make_axes_locatable
 
-        # get a map showing the reference or target image
+        # get background image (reference or target image)
         if backgroundIm not in ['tgt', 'ref']:
             raise ValueError('backgroundIm')
         backgroundIm = self.im2shift if backgroundIm == 'tgt' else self.imref
+
+        # make sure the output figure has a reasonable size, also if figsize is not given
+        if not figsize:
+            r, c = backgroundIm.shape[:2]
+            figsize = (8, r / c * 8) if r > c else (c / r * 8, 8)
+
+        # apply figsize multiplier
+        if figsize_multiplier:
+            if figsize_multiplier < 0:
+                raise ValueError(figsize_multiplier, 'The figure size multiplier must be a positive finite number.')
+            figsize = (figsize[0] * figsize_multiplier, figsize[1] * figsize_multiplier)
+
+        # get a map showing the background image
         fig, ax = backgroundIm.show_map(figsize=figsize,
                                         nodataVal=self.nodata[1],
                                         return_map=True,
                                         band=self.COREG_obj.shift.band4match)
 
-        # make sure the output figure has a reasonable size, also if figsize is not given
-        if not figsize:
-            w, h = fig.get_size_inches()
-            fig.set_size_inches(w * 1.6, h * 1.6)
-
         # set figure title
         dict_attr_title = dict(
             X_WIN_SIZE='size of the matching window in x-direction [pixels]',
@@ -689,8 +699,11 @@ class COREG_LOCAL(object):
             if not self.q:
                 warnings.warn(msg)
 
+        # remove white space around the figure
+        plt.subplots_adjust(top=.95, bottom=.05, right=.95, left=.05)
+
         if savefigPath:
-            fig.savefig(savefigPath, dpi=savefigDPI)
+            fig.savefig(savefigPath, dpi=savefigDPI, pad_inches=0.1, bbox_inches='tight')
 
         if return_map:
             return fig, ax

arosics/Tie_Point_Grid.py

@@ -492,6 +492,100 @@ class Tie_Point_Grid(object):
 
         return float(np.median(ssim_col))
 
+    def calc_overall_stats(self, include_outliers: bool = False) -> dict:
+        """Calculate statistics like RMSE, MSE, MAE, ... from the tie point grid.
+
+        Full list of returned statistics:
+
+        - N_TP:                 number of tie points
+        - N_VALID_TP:           number of valid tie points
+        - N_INVALID_TP:         number of invalid tie points (false-positives)
+        - PERC_VALID_TP:        percentage of valid tie points
+        - RMSE_M:               root mean squared error of absolute shift vector length in map units
+        - RMSE_X_M:             root mean squared error of shift vector length in x-direction in map units
+        - RMSE_Y_M:             root mean squared error of shift vector length in y-direction in map units
+        - RMSE_X_PX:            root mean squared error of shift vector length in x-direction in pixel units
+        - RMSE_Y_PX:            root mean squared error of shift vector length in y-direction in pixel units
+        - MSE_M:                mean squared error of absolute shift vector length in map units
+        - MSE_X_M:              mean squared error of shift vector length in x-direction in map units
+        - MSE_Y_M:              mean squared error of shift vector length in y-direction in map units
+        - MSE_X_PX:             mean squared error of shift vector length in x-direction in pixel units
+        - MSE_Y_PX:             mean squared error of shift vector length in y-direction in pixel units
+        - MAE_M:                mean absolute error of absolute shift vector length in map units
+        - MAE_X_M:              mean absolute error of shift vector length in x-direction in map units
+        - MAE_Y_M:              mean absolute error of shift vector length in y-direction in map units
+        - MAE_X_PX:             mean absolute error of shift vector length in x-direction in pixel units
+        - MAE_Y_PX:             mean absolute error of shift vector length in y-direction in pixel units
+        - MEAN_ANGLE:           mean direction of the shift vectors in degrees from north
+        - MEDIAN_ANGLE:         median direction of the shift vectors in degrees from north
+        - MEAN_SSIM_BEFORE:     mean structural similatity index within each matching window before co-registration
+        - MEDIAN_SSIM_BEFORE:   median structural similatity index within each matching window before co-registration
+        - MEAN_SSIM_AFTER:      mean structural similatity index within each matching window after co-registration
+        - MEDIAN_RELIABILITY:   median tie point reliability in percent
+        - MEDIAN_SSIM_AFTER:    median structural similatity index within each matching window after co-registration
+        - MEAN_RELIABILITY:     mean tie point reliability in percent
+
+        :param include_outliers:    whether to include tie points that have been marked as false-positives (if present)
+        """
+        if self.CoRegPoints_table.empty:
+            raise RuntimeError('Cannot compute overall statistics because no tie points were found at all.')
+
+        tbl = self.CoRegPoints_table
+
+        n_tiepoints = sum(tbl['ABS_SHIFT'] != self.outFillVal)
+        n_outliers = sum(tbl['OUTLIER'] == 1)
+
+        tbl = tbl if include_outliers else tbl[tbl['OUTLIER'] == 0].copy() if 'OUTLIER' in tbl.columns else tbl
+        tbl = tbl.copy().replace(self.outFillVal, np.nan)
+
+        def RMSE(vals):
+            vals_sq = vals ** 2
+            return np.sqrt(sum(vals_sq) / len(vals_sq))
+
+        def MSE(vals):
+            vals_sq = vals ** 2
+            return sum(vals_sq) / len(vals_sq)
+
+        def MAE(vals):
+            vals_abs = np.abs(vals)
+            return sum(vals_abs) / len(vals_abs)
+
+        abs_shift, x_shift_m, y_shift_m, x_shift_px, y_shift_px, angle, ssim_before, ssim_after, reliability = \
+            [tbl[k].dropna().values for k in ['ABS_SHIFT', 'X_SHIFT_M', 'Y_SHIFT_M', 'X_SHIFT_PX', 'Y_SHIFT_PX',
+                                              'ANGLE', 'SSIM_BEFORE', 'SSIM_AFTER', 'RELIABILITY']]
+
+        stats = dict(
+            N_TP=n_tiepoints,
+            N_VALID_TP=len(abs_shift),
+            N_INVALID_TP=n_outliers,
+            PERC_VALID_TP=(n_tiepoints - n_outliers) / n_tiepoints * 100,
+            RMSE_M=RMSE(abs_shift),
+            RMSE_X_M=RMSE(x_shift_m),
+            RMSE_Y_M=RMSE(y_shift_m),
+            RMSE_X_PX=RMSE(x_shift_px),
+            RMSE_Y_PX=RMSE(y_shift_px),
+            MSE_M=MSE(abs_shift),
+            MSE_X_M=MSE(x_shift_m),
+            MSE_Y_M=MSE(y_shift_m),
+            MSE_X_PX=MSE(x_shift_px),
+            MSE_Y_PX=MSE(y_shift_px),
+            MAE_M=MAE(abs_shift),
+            MAE_X_M=MAE(x_shift_m),
+            MAE_Y_M=MAE(y_shift_m),
+            MAE_X_PX=MAE(x_shift_px),
+            MAE_Y_PX=MAE(y_shift_px),
+            MEAN_ANGLE=np.mean(angle),
+            MEDIAN_ANGLE=np.median(angle),
+            MEAN_SSIM_BEFORE=np.mean(ssim_before),
+            MEDIAN_SSIM_BEFORE=np.median(ssim_before),
+            MEAN_SSIM_AFTER=np.mean(ssim_after),
+            MEDIAN_SSIM_AFTER=np.median(ssim_after),
+            MEAN_RELIABILITY=np.mean(reliability),
+            MEDIAN_RELIABILITY=np.median(reliability)
+        )
+
+        return stats
+
     def plot_shift_distribution(self,
                                 include_outliers: bool = True,
                                 unit: str = 'm',
@@ -500,18 +594,26 @@ class Tie_Point_Grid(object):
                                 xlim: list = None,
                                 ylim: list = None,
                                 fontsize: int = 12,
-                                title: str = 'shift distribution'
+                                title: str = 'shift distribution',
+                                savefigPath: str = '',
+                                savefigDPI: int = 96,
+                                showFig: bool = True,
+                                return_fig: bool = False
                                 ) -> tuple:
         """Create a 2D scatterplot containing the distribution of calculated X/Y-shifts.
 
         :param include_outliers:    whether to include tie points that have been marked as false-positives
         :param unit:                'm' for meters or 'px' for pixels (default: 'm')
-        :param interactive:         interactive mode uses plotly for visualization
+        :param interactive:         whether to use interactive mode (uses plotly for visualization)
         :param figsize:             (xdim, ydim)
         :param xlim:                [xmin, xmax]
         :param ylim:                [ymin, ymax]
         :param fontsize:            size of all used fonts
         :param title:               the title to be plotted above the figure
+        :param savefigPath:         path where to save the figure
+        :param savefigDPI:          DPI resolution of the output figure when saved to disk
+        :param showFig:             whether to show or to hide the figure
+        :param return_fig:          whether to return the figure and axis objects
         """
         from matplotlib import pyplot as plt
 
@@ -598,9 +700,19 @@ class Tie_Point_Grid(object):
             leg = plt.legend(reversed(handles), reversed(labels), fontsize=fontsize, loc='upper right', scatterpoints=3)
             leg.get_frame().set_edgecolor('black')
 
-            plt.show()
+            # remove white space around the figure
+            plt.subplots_adjust(top=.94, bottom=.06, right=.96, left=.09)
 
-            return fig, ax
+            if savefigPath:
+                fig.savefig(savefigPath, dpi=savefigDPI, pad_inches=0.3, bbox_inches='tight')
+
+            if return_fig:
+                return fig, ax
+
+            if showFig and not self.q:
+                plt.show(block=True)
+            else:
+                plt.close(fig)
 
     def dump_CoRegPoints_table(self, path_out=None):
         if self.CoRegPoints_table.empty:
@@ -685,7 +797,7 @@ class Tie_Point_Grid(object):
                           skip_nodata: bool = True,
                           skip_nodata_col: str = 'ABS_SHIFT'
                           ) -> None:
-        """Write the calculated tie points grid to a point shapefile (e.g., for visualization by a GIS software).
+        """Write the calculated tie point grid to a point shapefile (e.g., for visualization by a GIS software).
 
         NOTE: The shapefile uses Tie_Point_Grid.CoRegPoints_table as attribute table.
 

arosics/version.py

@@ -24,5 +24,5 @@
 # limitations under the License.
 
 
-__version__ = '1.6.0'
-__versionalias__ = '2021-09-27_01'
+__version__ = '1.7.0'
+__versionalias__ = '2021-09-30_01'

tests/test_tie_point_grid.py

@@ -80,6 +80,16 @@ class Test_Tie_Point_Grid(unittest.TestCase):
         self.TPG.calc_overall_ssim(include_outliers=False, after_correction=True)
         self.TPG.calc_overall_ssim(include_outliers=True, after_correction=False)
 
+    def test_calc_overall_stats(self):
+        stats_noOL = self.TPG.calc_overall_stats(include_outliers=False)
+        stats_OL = self.TPG.calc_overall_stats(include_outliers=True)
+
+        self.assertTrue(stats_noOL)
+        self.assertTrue(stats_OL)
+        self.assertIsInstance(stats_noOL, dict)
+        self.assertIsInstance(stats_OL, dict)
+        self.assertNotEqual(stats_noOL, stats_OL)
+
     def test_plot_shift_distribution(self):
         with warnings.catch_warnings():
             warnings.filterwarnings(