Coverage for arosics/CoReg.py: 86%
735 statements
« prev ^ index » next coverage.py v7.3.2, created at 2024-04-03 14:59 +0000
1# -*- coding: utf-8 -*-
3# AROSICS - Automated and Robust Open-Source Image Co-Registration Software
4#
5# Copyright (C) 2017-2024
6# - Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de)
7# - Helmholtz Centre Potsdam - GFZ German Research Centre for Geosciences Potsdam,
8# Germany (https://www.gfz-potsdam.de/)
9#
10# This software was developed within the context of the GeoMultiSens project funded
11# by the German Federal Ministry of Education and Research
12# (project grant code: 01 IS 14 010 A-C).
13#
14# Licensed under the Apache License, Version 2.0 (the "License");
15# you may not use this file except in compliance with the License.
16# You may obtain a copy of the License at
17#
18# https://www.apache.org/licenses/LICENSE-2.0
19#
20# Unless required by applicable law or agreed to in writing, software
21# distributed under the License is distributed on an "AS IS" BASIS,
22# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
23# See the License for the specific language governing permissions and
24# limitations under the License.
26import os
27import time
28import warnings
29from copy import copy
30from typing import Iterable, Union, Tuple, List, Optional # noqa F401
31from multiprocessing import cpu_count
33# custom
34from osgeo import gdal # noqa
35import numpy as np
36from scipy.fft import fft2, ifft2, fftshift
38from packaging.version import parse as parse_version
39try:
40 import pyfftw
41 # pyfftw>=0.13.0 is currently not used due to https://github.com/pyFFTW/pyFFTW/issues/294
42 if parse_version(pyfftw.__version__) >= parse_version('0.13.0'):
43 pyfftw = None
44except ImportError:
45 pyfftw = None
46from shapely.geometry import Point, Polygon
48# internal modules
49from .DeShifter import DESHIFTER, _dict_rspAlg_rsp_Int
50from . import geometry as GEO
51from . import plotting as PLT
53from geoarray import GeoArray
54from py_tools_ds.convenience.object_oriented import alias_property
55from py_tools_ds.geo.coord_calc import get_corner_coordinates
56from py_tools_ds.geo.vector.topology import get_overlap_polygon, get_smallest_boxImYX_that_contains_boxMapYX
57from py_tools_ds.geo.projection import prj_equal
58from py_tools_ds.geo.vector.geometry import boxObj, round_shapelyPoly_coords
59from py_tools_ds.geo.coord_grid import move_shapelyPoly_to_image_grid, is_coord_grid_equal
60from py_tools_ds.geo.coord_trafo import reproject_shapelyGeometry, mapXY2imXY, imXY2mapXY
61from py_tools_ds.geo.raster.reproject import warp_ndarray
62from py_tools_ds.geo.map_info import geotransform2mapinfo
63from py_tools_ds.io.vector.writer import write_shp
65__author__ = 'Daniel Scheffler'
68class GeoArray_CoReg(GeoArray):
69 def __init__(self,
70 CoReg_params: dict,
71 imID: str
72 ) -> None:
73 assert imID in ['ref', 'shift']
75 # run GeoArray init
76 path_or_geoArr = CoReg_params['im_ref'] if imID == 'ref' else CoReg_params['im_tgt']
77 nodata = CoReg_params['nodata'][0 if imID == 'ref' else 1]
78 progress = CoReg_params['progress']
79 q = CoReg_params['q'] if not CoReg_params['v'] else False
81 super(GeoArray_CoReg, self).__init__(path_or_geoArr, nodata=nodata, progress=progress, q=q)
83 self.imID = imID
84 self.imName = 'reference image' if imID == 'ref' else 'image to be shifted'
85 self.v = CoReg_params['v']
87 assert isinstance(self, GeoArray), \
88 'Something went wrong with the creation of GeoArray instance for the %s. The created ' \
89 'instance does not seem to belong to the GeoArray class. If you are working in Jupyter Notebook, reset ' \
90 'the kernel and try again.' % self.imName
92 # set title to be used in plots
93 self.title = os.path.basename(self.filePath) if self.filePath else self.imName
95 # validate params
96 # assert self.prj, 'The %s has no projection.' % self.imName # TODO
97 # assert not re.search('LOCAL_CS', self.prj), 'The %s is not georeferenced.' % self.imName # TODO
98 assert self.gt, 'The %s has no map information.' % self.imName
100 # set band4match
101 self.band4match = (CoReg_params['r_b4match'] if imID == 'ref' else CoReg_params['s_b4match']) - 1
102 assert self.bands >= self.band4match + 1 >= 1, \
103 "The %s has %s %s. So its band number to match must be %s%s. Got %s." \
104 % (self.imName, self.bands, 'bands' if self.bands > 1 else
105 'band', 'between 1 and ' if self.bands > 1 else '', self.bands, self.band4match)
107 # compute nodata mask and validate that it is not completely filled with nodata
108 self.calc_mask_nodata(fromBand=self.band4match) # this avoids that all bands have to be read
110 if True not in self.mask_nodata[:]:
111 raise RuntimeError(f'The {self.imName} passed to AROSICS only contains nodata values.')
113 # set footprint_poly
114 given_footprint_poly = CoReg_params['footprint_poly_%s' % ('ref' if imID == 'ref' else 'tgt')]
115 given_corner_coord = CoReg_params['data_corners_%s' % ('ref' if imID == 'ref' else 'tgt')]
117 if given_footprint_poly:
118 self.footprint_poly = given_footprint_poly
119 elif given_corner_coord is not None:
120 self.footprint_poly = Polygon(given_corner_coord)
121 elif not CoReg_params['calc_corners']:
122 # use the image extent
123 self.footprint_poly = Polygon(get_corner_coordinates(gt=self.gt, cols=self.cols, rows=self.rows))
124 else:
125 # footprint_poly is calculated automatically by GeoArray
126 if not CoReg_params['q']:
127 print('Calculating footprint polygon and actual data corner coordinates for %s...' % self.imName)
129 with warnings.catch_warnings(record=True) as w:
130 _ = self.footprint_poly # execute getter
132 if len(w) > 0 and 'disjunct polygon(s) outside' in str(w[-1].message):
133 warnings.warn('The footprint of the %s contains multiple separate image parts. '
134 'AROSICS will only process the largest image part.' % self.imName)
135 # FIXME use a convex hull as footprint poly
137 # validate footprint poly
138 if not self.footprint_poly.is_valid:
139 self.footprint_poly = self.footprint_poly.buffer(0)
141 if not self.q:
142 print('Bounding box of calculated footprint for %s:\n\t%s' % (self.imName, self.footprint_poly.bounds))
144 # add bad data mask
145 given_mask = CoReg_params['mask_baddata_%s' % ('ref' if imID == 'ref' else 'tgt')]
146 if given_mask is not None:
147 self.mask_baddata = given_mask # runs GeoArray.mask_baddata.setter -> sets it to BadDataMask()
149 poly = alias_property('footprint_poly') # ensures that self.poly is updated if self.footprint_poly is updated
152class COREG(object):
153 """The COREG class detects and corrects global X/Y shifts between a target and reference image.
155 Geometric shifts are calculated at a specific (adjustable) image position. Correction performs a global shifting
156 in X- or Y direction.
158 See help(COREG) for documentation!
159 """
161 def __init__(self,
162 im_ref: Union[GeoArray, str],
163 im_tgt: Union[GeoArray, str],
164 path_out: str = None,
165 fmt_out: str = 'ENVI',
166 out_crea_options: list = None,
167 r_b4match: int = 1,
168 s_b4match: int = 1,
169 wp: Tuple[float, float] = (None, None),
170 ws: Tuple[int, int] = (256, 256),
171 max_iter: int = 5,
172 max_shift: int = 5,
173 align_grids: bool = False,
174 match_gsd: bool = False,
175 out_gsd: Tuple[float] = None,
176 target_xyGrid: List[List] = None,
177 resamp_alg_deshift: str = 'cubic',
178 resamp_alg_calc: str = 'cubic',
179 footprint_poly_ref: str = None,
180 footprint_poly_tgt: str = None,
181 data_corners_ref: list = None,
182 data_corners_tgt: list = None,
183 nodata: Tuple = (None, None),
184 calc_corners: bool = True,
185 binary_ws: bool = True,
186 mask_baddata_ref: Union[GeoArray, str] = None,
187 mask_baddata_tgt: Union[GeoArray, str] = None,
188 CPUs: int = None,
189 force_quadratic_win: bool = True,
190 progress: bool = True,
191 v: bool = False,
192 path_verbose_out: str = None,
193 q: bool = False,
194 ignore_errors: bool = False
195 ) -> None:
196 """Get an instance of the COREG class.
198 :param im_ref:
199 source path or GeoArray instance of reference image
200 (any GDAL compatible image format is supported)
202 :param im_tgt:
203 source path or GeoArray instance of the target image, i.e., the image to be shifted
204 (any GDAL compatible image format is supported)
206 :param path_out:
207 target path of the coregistered image
208 - if None (default), the method correct_shifts() does not write to disk
209 - if 'auto': /dir/of/im1/<im1>__shifted_to__<im0>.bsq
211 :param fmt_out:
212 raster file format for output file. ignored if path_out is None. can be any GDAL compatible raster file
213 format (e.g. 'ENVI', 'GTIFF'; default: ENVI). Refer to https://gdal.org/drivers/raster/index.html to get a
214 full list of supported formats.
216 :param out_crea_options:
217 GDAL creation options for the output image, e.g. ["QUALITY=80", "REVERSIBLE=YES", "WRITE_METADATA=YES"]
219 :param r_b4match:
220 band of reference image to be used for matching (starts with 1; default: 1)
222 :param s_b4match:
223 band of shift image to be used for matching (starts with 1; default: 1)
225 :param wp:
226 custom matching window position as (X, Y) map coordinate in the same projection as the reference image
227 (default: central position of image overlap)
229 :param ws:
230 custom matching window size [pixels] as (X, Y) tuple (default: (256,256))
232 :param max_iter:
233 maximum number of iterations for matching (default: 5)
235 :param max_shift:
236 maximum shift distance in reference image pixel units (default: 5 px)
238 :param align_grids:
239 True: align the input coordinate grid to the reference (does not affect the output pixel size as long as
240 input and output pixel sizes are compatible (5:30 or 10:30 but not 4:30)), default = False
242 - NOTE: If this is set to False, the mis-registration in the target image is corrected by updating its
243 geocoding information only, i.e., without performing any resampling which preserves the original
244 image values (except that a resampling is needed due to other reasons, e.g., if 'match_gsd',
245 'out_gsd' or 'target_xyGrid' are given).
247 :param match_gsd:
248 True: match the input pixel size to the reference pixel size,
249 default = False
251 :param out_gsd:
252 output pixel size (X/Y) in units of the reference coordinate system (default = pixel size of the target
253 image), given values are overridden by match_gsd=True
255 :param target_xyGrid:
256 a list with a target x-grid and a target y-grid like [[15,45], [15,45]]
257 This overrides 'out_gsd', 'align_grids' and 'match_gsd'.
259 :param resamp_alg_deshift:
260 the resampling algorithm to be used for shift correction (if neccessary)
261 - valid algorithms: nearest, bilinear, cubic, cubic_spline, lanczos, average, mode, max, min, med, q1, q3
262 - default: cubic
264 :param resamp_alg_calc:
265 the resampling algorithm to be used for all warping processes during calculatio of spatial shift
266 - valid algorithms: nearest, bilinear, cubic, cubic_spline, lanczos, average, mode, max, min, med, q1, q3
267 - default: cubic (highly recommended)
269 :param footprint_poly_ref:
270 footprint polygon of the reference image (WKT string or shapely.geometry.Polygon),
271 e.g. 'POLYGON ((299999 6000000, 299999 5890200, 409799 5890200, 409799 6000000, 299999 6000000))'
273 :param footprint_poly_tgt:
274 footprint polygon of the image to be shifted (WKT string or shapely.geometry.Polygon)
275 e.g. 'POLYGON ((299999 6000000, 299999 5890200, 409799 5890200, 409799 6000000, 299999 6000000))'
277 :param data_corners_ref:
278 map coordinates of data corners within reference image. ignored if footprint_poly_ref is given.
280 :param data_corners_tgt:
281 map coordinates of data corners within image to be shifted. ignored if footprint_poly_tgt is given.
283 :param nodata:
284 no-data values for reference image and image to be shifted. The default is (None, None) which indicates
285 that there is no specific no-data value for both of the input images.
287 :param calc_corners:
288 calculate true positions of the dataset corners in order to get a useful matching window position within
289 the actual image overlap (default: True; deactivated if '-cor0' and '-cor1' are given)
291 :param binary_ws:
292 use binary X/Y dimensions for the matching window (default: True)
294 :param mask_baddata_ref:
295 path to a 2D boolean mask file (or an instance of GeoArray) for the reference image where all bad data
296 pixels (e.g. clouds) are marked with True and the remaining pixels with False. Must have the same
297 geographic extent and projection like 'im_ref'. The mask is used to check if the chosen matching window
298 position is valid in the sense of useful data. Otherwise, this window position is rejected.
300 :param mask_baddata_tgt:
301 path to a 2D boolean mask file (or an instance of GeoArray) for the image to be shifted where all bad data
302 pixels (e.g. clouds) are marked with True and the remaining pixels with False. Must have the same
303 geographic extent and projection like 'im_ref'. The mask is used to check if the chosen matching window
304 position is valid in the sense of useful data. Otherwise, this window position is rejected.
306 :param CPUs:
307 number of CPUs to use during pixel grid equalization (default: None, which means 'all CPUs available')
309 :param force_quadratic_win:
310 force a quadratic matching window (default: True)
312 :param progress:
313 show progress bars (default: True)
315 :param v:
316 verbose mode (default: False)
318 :param path_verbose_out:
319 an optional output directory for intermediate results
320 (if not given, no intermediate results are written to disk)
322 :param q:
323 quiet mode (default: False)
325 :param ignore_errors:
326 Useful for batch processing. (default: False)
327 In case of error COREG.success == False and COREG.x_shift_px/COREG.y_shift_px is None
328 """
329 self.params = dict([x for x in locals().items() if x[0] != "self"])
331 # input validation
332 if gdal.GetDriverByName(fmt_out) is None:
333 raise ValueError(fmt_out, "'%s' is not a supported GDAL driver." % fmt_out)
335 if match_gsd and out_gsd:
336 warnings.warn("'-out_gsd' is ignored because '-match_gsd' is set.\n")
338 if out_gsd and (not isinstance(out_gsd, list) or len(out_gsd) != 2):
339 raise ValueError(out_gsd, 'out_gsd must be a list with two values.')
341 if data_corners_ref and not isinstance(data_corners_ref[0], list):
342 # group if not [[x,y],[x,y]..] but [x,y,x,y,]
343 data_corners_ref = [data_corners_ref[i:i + 2]
344 for i in range(0, len(data_corners_ref), 2)]
346 if data_corners_tgt and not isinstance(data_corners_tgt[0], list):
347 # group if not [[x,y],[x,y]..]
348 data_corners_tgt = [data_corners_tgt[i:i + 2]
349 for i in range(0, len(data_corners_tgt), 2)]
351 if nodata and (not isinstance(nodata, Iterable) or len(nodata) != 2):
352 raise ValueError(nodata, "'nodata' must be an iterable with two values. "
353 "Got %s with length %s." % (type(nodata), len(nodata)))
355 for rspAlg in [resamp_alg_deshift, resamp_alg_calc]:
356 if rspAlg not in _dict_rspAlg_rsp_Int.keys():
357 raise ValueError("'%s' is not a supported resampling algorithm." % rspAlg)
359 if resamp_alg_calc in ['average', 5] and (v or not q):
360 warnings.warn("The resampling algorithm 'average' causes sinus-shaped patterns in fft images that will "
361 "affect the precision of the calculated spatial shifts! It is highly recommended to "
362 "choose another resampling algorithm.")
364 self.path_out = path_out # updated by self.set_outpathes
365 self.fmt_out = fmt_out
366 self.out_creaOpt = out_crea_options
367 self.win_pos_XY = wp # updated by self.get_opt_winpos_winsize()
368 self.win_size_XY = ws # updated by self.get_opt_winpos_winsize()
369 self.max_iter = max_iter
370 self.max_shift = max_shift
371 self.align_grids = align_grids
372 self.match_gsd = match_gsd
373 self.out_gsd = out_gsd
374 self.target_xyGrid = target_xyGrid
375 self.rspAlg_DS = resamp_alg_deshift \
376 if isinstance(resamp_alg_deshift, str) else _dict_rspAlg_rsp_Int[resamp_alg_deshift]
377 self.rspAlg_calc = resamp_alg_calc \
378 if isinstance(resamp_alg_calc, str) else _dict_rspAlg_rsp_Int[resamp_alg_calc]
379 self.calc_corners = calc_corners
380 self.CPUs = CPUs if CPUs and CPUs <= cpu_count() else cpu_count()
381 self.bin_ws = binary_ws
382 self.force_quadratic_win = force_quadratic_win
383 self.v = v
384 self.path_verbose_out = path_verbose_out
385 self.q = q if not v else False # overridden by v
386 self.progress = progress if not q else False # overridden by q
388 self.ignErr = ignore_errors
389 self.max_win_sz_changes = 3 # TODO: änderung der window size, falls nach max_iter kein valider match gefunden
390 self.ref: Optional[GeoArray_CoReg] = None # set by self.get_image_params
391 self.shift: Optional[GeoArray_CoReg] = None # set by self.get_image_params
392 self.matchBox: Optional[boxObj] = None # set by self.get_clip_window_properties()
393 self.otherBox: Optional[boxObj] = None # set by self.get_clip_window_properties()
394 self.matchWin: Optional[GeoArray] = None # set by self._get_image_windows_to_match()
395 self.otherWin: Optional[GeoArray] = None # set by self._get_image_windows_to_match()
396 self.overlap_poly = None # set by self._get_overlap_properties()
397 self.overlap_percentage = None # set by self._get_overlap_properties()
398 self.overlap_area = None # set by self._get_overlap_properties()
399 self.imfft_xgsd = None # set by self.get_clip_window_properties()
400 self.imfft_ygsd = None # set by self.get_clip_window_properties()
401 self.fftw_works = None # set by self._calc_shifted_cross_power_spectrum()
402 self.fftw_win_size_YX = None # set by calc_shifted_cross_power_spectrum()
404 self.x_shift_px = None # always in shift image units (image coords) # set by calculate_spatial_shifts()
405 self.y_shift_px = None # always in shift image units (image coords) # set by calculate_spatial_shifts()
406 self.x_shift_map = None # set by self.get_updated_map_info()
407 self.y_shift_map = None # set by self.get_updated_map_info()
408 self.vec_length_map = None
409 self.vec_angle_deg = None
410 self.updated_map_info = None # set by self.get_updated_map_info()
411 self.ssim_orig = None # set by self._validate_ssim_improvement()
412 self.ssim_deshifted = None # set by self._validate_ssim_improvement()
413 self._ssim_improved = None # private attribute to be filled by self.ssim_improved
414 self.shift_reliability = None # set by self.calculate_spatial_shifts()
416 self.tracked_errors = [] # expanded each time an error occurs
417 self.success = None # default
418 self.deshift_results = None # set by self.correct_shifts()
420 self._check_and_handle_metaRotation()
421 self._get_image_params()
422 self._set_outpathes(im_ref, im_tgt)
423 self.grid2use = 'ref' if self.shift.xgsd <= self.ref.xgsd else 'shift'
424 if self.v:
425 print('resolutions: ', self.ref.xgsd, self.shift.xgsd)
427 self._get_overlap_properties()
429 if self.v and self.path_verbose_out:
430 write_shp(os.path.join(self.path_verbose_out, 'poly_imref.shp'), self.ref.poly, self.ref.prj)
431 write_shp(os.path.join(self.path_verbose_out, 'poly_im2shift.shp'), self.shift.poly, self.shift.prj)
432 write_shp(os.path.join(self.path_verbose_out, 'overlap_poly.shp'), self.overlap_poly, self.ref.prj)
434 # FIXME: transform_mapPt1_to_mapPt2(im2shift_center_map, ds_imref.GetProjection(), ds_im2shift.GetProjection())
435 # FIXME später basteln für den fall, dass projektionen nicht gleich sind
437 # get_clip_window_properties
438 self._get_opt_winpos_winsize()
439 if not self.q:
440 print('Matching window position (X,Y): %s/%s' % (self.win_pos_XY[0], self.win_pos_XY[1]))
441 self._get_clip_window_properties() # sets self.matchBox, self.otherBox and much more
443 if self.v and self.path_verbose_out and self.matchBox.mapPoly and self.success is not False:
444 write_shp(os.path.join(self.path_verbose_out, 'poly_matchWin.shp'),
445 self.matchBox.mapPoly, self.matchBox.prj)
447 self.success = False if self.success is False or not self.matchBox.boxMapYX else None
449 self._coreg_info = None # private attribute to be filled by self.coreg_info property
451 def _handle_error(self, error, warn=False, warnMsg=None):
452 """Append the given error to self.tracked_errors.
454 This sets self.success to False and raises the error in case self.ignore_errors = True.
456 :param error: instance of an error
457 :param warn: whether to give a warning in case error would be ignored otherwise
458 :param warnMsg: a custom message for the warning
459 :return:
460 """
461 warn = warn or warnMsg is not None or self.v
463 self.tracked_errors.append(error)
464 self.success = False
466 if self.ignErr and warn:
467 warnMsg = repr(error) if not warnMsg else warnMsg
468 print('\nWARNING: ' + warnMsg)
470 if not self.ignErr:
471 raise error
473 def _set_outpathes(self,
474 im_ref: Union[(GeoArray, str)],
475 im_tgt: Union[(GeoArray, str)]
476 ) -> None:
477 def get_baseN(path):
478 return os.path.splitext(os.path.basename(path))[0]
480 # get input paths
481 def get_input_path(im):
482 path = im.filePath if isinstance(im, GeoArray) else im
484 if isinstance(im, GeoArray) and im.filePath is None and self.path_out == 'auto':
485 raise ValueError(self.path_out, "The output path must be explicitly set in case the input "
486 "reference or target image is in-memory (without a reference to a "
487 "physical file on disk). Received path_out='%s'." % self.path_out)
489 return path
491 path_im_ref = get_input_path(im_ref)
492 path_im_tgt = get_input_path(im_tgt)
494 if self.path_out: # this also applies to self.path_out='auto'
496 if self.path_out == 'auto':
497 dir_out, fName_out = os.path.dirname(path_im_tgt), ''
498 else:
499 dir_out, fName_out = os.path.split(self.path_out)
501 if dir_out and fName_out:
502 # a valid output path is given => do nothing
503 pass
505 else:
506 # automatically create an output directory and filename if not given
507 if not dir_out:
508 if not path_im_ref:
509 dir_out = os.path.abspath(os.path.curdir)
510 else:
511 dir_out = os.path.dirname(path_im_ref)
513 if not fName_out:
514 ext = 'bsq' if self.fmt_out == 'ENVI' else \
515 gdal.GetDriverByName(self.fmt_out).GetMetadataItem(gdal.DMD_EXTENSION)
516 fName_out = fName_out if fName_out not in ['.', ''] else \
517 '%s__shifted_to__%s' % (get_baseN(path_im_tgt), get_baseN(path_im_ref))
518 fName_out = fName_out + '.%s' % ext if ext else fName_out
520 self.path_out = os.path.abspath(os.path.join(dir_out, fName_out))
522 assert ' ' not in self.path_out, \
523 "The path of the output image contains whitespaces. This is not supported by GDAL."
524 else:
525 # this only happens if COREG is not instanced from within Python and self.path_out is explicitly set to None
526 # => DESHIFTER will return an array
527 pass
529 if self.v:
530 if self.path_verbose_out:
531 dir_out, dirname_out = os.path.split(self.path_verbose_out)
533 if not dir_out:
534 if self.path_out:
535 self.path_verbose_out = os.path.dirname(self.path_out)
536 else:
537 self.path_verbose_out = \
538 os.path.abspath(os.path.join(os.path.curdir, 'CoReg_verboseOut__%s__shifted_to__%s'
539 % (get_baseN(path_im_tgt), get_baseN(path_im_ref))))
540 elif dirname_out and not dir_out:
541 self.path_verbose_out = os.path.abspath(os.path.join(os.path.curdir, dirname_out))
543 assert ' ' not in self.path_verbose_out, \
544 "'path_verbose_out' contains whitespaces. This is not supported by GDAL."
546 else:
547 self.path_verbose_out = None
549 if self.path_verbose_out and not os.path.isdir(self.path_verbose_out):
550 os.makedirs(self.path_verbose_out)
552 def _get_image_params(self) -> None:
553 self.ref = GeoArray_CoReg(self.params, 'ref')
554 self.shift = GeoArray_CoReg(self.params, 'shift')
556 if not prj_equal(self.ref.prj, self.shift.prj):
557 from pyproj import CRS
559 crs_ref = CRS.from_user_input(self.ref.prj)
560 crs_shift = CRS.from_user_input(self.shift.prj)
562 name_ref, name_shift = \
563 (crs_ref.name, crs_shift.name) if not crs_ref.name == crs_shift.name else (crs_ref.srs, crs_shift.srs)
565 raise RuntimeError(
566 'Input projections are not equal. Different projections are currently not supported. '
567 'Got %s vs. %s.' % (name_ref, name_shift))
569 def _get_overlap_properties(self) -> None:
570 with warnings.catch_warnings():
571 # already warned in GeoArray_CoReg.__init__()
572 warnings.filterwarnings("ignore", category=RuntimeWarning, message=".*disjunct.*")
574 overlap_tmp = get_overlap_polygon(self.ref.poly, self.shift.poly, self.v)
576 self.overlap_poly = overlap_tmp['overlap poly'] # has to be in reference projection
577 self.overlap_percentage = overlap_tmp['overlap percentage']
578 self.overlap_area = overlap_tmp['overlap area']
580 assert self.overlap_poly, 'The input images have no spatial overlap.'
582 # overlap are must at least cover 16*16 pixels
583 px_area = self.ref.xgsd * self.ref.ygsd if self.grid2use == 'ref' else self.shift.xgsd * self.shift.ygsd
584 px_covered = self.overlap_area / px_area
585 assert px_covered > 16 * 16, \
586 'Overlap area covers only %s pixels. At least 16*16 pixels are needed.' % px_covered
588 def _check_and_handle_metaRotation(self):
589 """Check if the provided input data have a metadata rotation and if yes, correct it.
591 In case there is a rotation, the GDAL GeoTransform is not 0 at positions 2 or 4. So far, AROSICS does not
592 handle such rotations, so the resampling is needed to make things work.
593 """
594 gA_ref = GeoArray(self.params['im_ref'])
595 gA_tgt = GeoArray(self.params['im_tgt'])
597 msg = 'The %s image needs to be resampled because it has a row/column rotation in ' \
598 'its map info which is not handled by AROSICS.'
600 if GEO.has_metaRotation(gA_ref):
601 warnings.warn(msg % 'reference')
602 self.params['im_ref'] = GEO.remove_metaRotation(gA_ref)
604 if GEO.has_metaRotation(gA_tgt):
605 warnings.warn(msg % 'target')
606 self.params['im_tgt'] = GEO.remove_metaRotation(gA_tgt)
608 @property
609 def are_pixGrids_equal(self):
610 return prj_equal(self.ref.prj, self.shift.prj) and \
611 is_coord_grid_equal(self.ref.gt, *self.shift.xygrid_specs, tolerance=1e-8)
613 def equalize_pixGrids(self) -> None:
614 """Equalize image grids and projections of reference and target image (align target to reference).
616 NOTE: This method is only called by COREG_LOCAL to speed up processing during detection of displacements.
617 """
618 if not self.are_pixGrids_equal or GEO.has_metaRotation(self.ref) or GEO.has_metaRotation(self.shift):
619 if not self.q:
620 print("Equalizing pixel grids and projections of reference and target image...")
622 def equalize(gA_from: GeoArray, gA_to: GeoArray) -> GeoArray:
623 if gA_from.bands > 1:
624 gA_from = gA_from.get_subset(zslice=slice(gA_from.band4match, gA_from.band4match + 1))
625 gA_from.reproject_to_new_grid(prototype=gA_to, CPUs=self.CPUs)
626 gA_from.band4match = 0 # after resampling there is only one band in the GeoArray
628 return gA_from
630 if self.grid2use == 'ref':
631 # resample target to reference image
632 self.shift = equalize(gA_from=self.shift, gA_to=self.ref)
634 else:
635 # resample reference to target image
636 self.ref = equalize(gA_from=self.ref, gA_to=self.shift)
638 # self.ref.gt = (self.ref.gt[0], 1, self.ref.gt[2], self.ref.gt[3], self.ref.gt[4], -1)
639 # self.shift.gt = (self.shift.gt[0], 1, self.shift.gt[2], self.shift.gt[3], self.shift.gt[4], -1)
641 def show_image_footprints(self):
642 """Show a web map containing the calculated footprints and overlap area of the input images.
644 NOTE: This method is intended to be called from Jupyter Notebook.
645 """
646 # TODO different colors for polygons
647 assert self.overlap_poly, 'Please calculate the overlap polygon first.'
649 import folium
650 import geojson
652 refPoly = reproject_shapelyGeometry(self.ref.poly, self.ref.prj, 4326)
653 shiftPoly = reproject_shapelyGeometry(self.shift.poly, self.shift.prj, 4326)
654 overlapPoly = reproject_shapelyGeometry(self.overlap_poly, self.shift.prj, 4326)
655 matchBoxPoly = reproject_shapelyGeometry(self.matchBox.mapPoly, self.shift.prj, 4326)
657 m = folium.Map(location=tuple(np.array(overlapPoly.centroid.coords.xy).flatten())[::-1])
658 for poly in [refPoly, shiftPoly, overlapPoly, matchBoxPoly]:
659 gjs = geojson.Feature(geometry=poly, properties={})
660 folium.GeoJson(gjs).add_to(m)
661 return m
663 def show_matchWin(self,
664 figsize: tuple = (15, 15),
665 interactive: bool = True,
666 after_correction: bool = None,
667 pmin=2,
668 pmax=98):
669 """Show the image content within the matching window.
671 :param figsize: figure size
672 :param interactive: whether to return an interactive figure based on 'holoviews' library
673 :param after_correction: True/False: show the image content AFTER shift correction or before
674 None: show both states - before and after correction (default)
675 :param pmin: percentage to be used for excluding the darkest pixels from stretching (default: 2)
676 :param pmax: percentage to be used for excluding the brightest pixels from stretching
677 (default: 98)
678 :return:
679 """
680 if not self.success and after_correction in [True, None]:
681 warnings.warn('It is only possible to show the matching window before correction of spatial displacements '
682 'because no valid displacement has been calculated yet.')
683 after_correction = False
685 if interactive:
686 # use Holoviews
687 try:
688 import holoviews as hv
689 except ImportError:
690 hv = None
691 if not hv:
692 raise ImportError(
693 "This method requires the library 'holoviews'. It can be installed for Conda with "
694 "the shell command 'conda install -c conda-forge holoviews bokeh'.")
696 hv.notebook_extension('matplotlib')
697 hv.Store.add_style_opts(hv.Image, ['vmin', 'vmax'])
699 # hv.Store.option_setters.options().Image = hv.Options('style', cmap='gnuplot2')
700 # hv.Store.add_style_opts(hv.Image, ['cmap'])
701 # renderer = hv.Store.renderers['matplotlib'].instance(fig='svg', holomap='gif')
702 # RasterPlot = renderer.plotting_class(hv.Image)
703 # RasterPlot.cmap = 'gray'
704 otherWin_corr = self._get_deshifted_otherWin() if after_correction in [True, None] else None
705 xmin, xmax, ymin, ymax = self.matchBox.boundsMap
707 def get_hv_image(geoArr: GeoArray):
708 from skimage.exposure import rescale_intensity # import here to avoid static TLS ImportError
710 arr_masked = np.ma.masked_equal(geoArr[:], geoArr.nodata)
711 vmin = np.nanpercentile(arr_masked.compressed(), pmin)
712 vmax = np.nanpercentile(arr_masked.compressed(), pmax)
713 arr2plot = rescale_intensity(arr_masked, in_range=(vmin, vmax), out_range='int8')
715 return (
716 hv.Image(
717 arr2plot,
718 bounds=(xmin, ymin, xmax, ymax)
719 ).options(
720 cmap='gray',
721 vmin=vmin,
722 vmax=vmax,
723 interpolation='none',
724 fig_inches=figsize,
725 show_grid=True
726 )
727 )
729 hvIm_matchWin = get_hv_image(self.matchWin)
730 hvIm_otherWin_orig = get_hv_image(self.otherWin)
731 hvIm_otherWin_corr = get_hv_image(otherWin_corr) if after_correction in [True, None] else None
733 if after_correction is None:
734 # view both states
735 print('Matching window before and after correction (left and right): ')
737 # get layouts (docs on options: https://holoviews.org/user_guide)
738 layout_before = (
739 hvIm_matchWin.options(title='BEFORE co-registration') +
740 hvIm_matchWin.options(title='AFTER co-registration')
741 ).options(
742 fig_inches=figsize
743 )
744 layout_after = (
745 hvIm_otherWin_orig.options(title='BEFORE co-registration') +
746 hvIm_otherWin_corr.options(title='AFTER co-registration')
747 ).options(
748 fig_inches=figsize
749 )
751 # plot!
752 imgs = {1: layout_before, 2: layout_after}
753 hmap = hv.HoloMap(imgs, kdims=['image']).collate().cols(2)
755 else:
756 # view state before or after correction
757 imgs = {1: hvIm_matchWin, 2: hvIm_otherWin_corr if after_correction else hvIm_otherWin_orig}
758 hmap = hv.HoloMap(imgs, kdims=['image'])
760 # Construct a HoloMap by evaluating the function over all the keys
761 # hmap = hv.HoloMap(imgs_corr, kdims=['image']) + hv.HoloMap(imgs_corr, kdims=['image'])
763 # Construct a HoloMap by defining the sampling on the Dimension
764 # dmap = hv.DynamicMap(image_slice, kdims=[hv.Dimension('z_axis', values=keys)])
766 return hmap
768 else:
769 # TODO add titles
770 # TODO handle after_correction=None here
771 self.matchWin.show(figsize=figsize)
772 if after_correction:
773 self._get_deshifted_otherWin().show(figsize=figsize, pmin=pmin, pmax=pmax)
774 else:
775 self.otherWin.show(figsize=figsize, pmin=pmin, pmax=pmax)
777 def show_cross_power_spectrum(self, interactive: bool = False) -> None:
778 """Show a 3D surface of the cross power spectrum.
780 NOTE: The cross power spectrum is the result from phase correlating the reference and target
781 image within the matching window.
783 :param interactive: whether to return an interactice 3D surface plot based on 'plotly' library
784 :return:
785 """
786 if interactive:
787 # create plotly 3D surface
789 # import plotly.plotly as py # online mode -> every plot is uploaded into online plotly account
790 from plotly.offline import iplot, init_notebook_mode
791 import plotly.graph_objs as go
793 init_notebook_mode(connected=True)
795 z_data = self._calc_shifted_cross_power_spectrum()
796 data = [go.Surface(z=z_data)]
797 layout = go.Layout(
798 title='cross power spectrum',
799 autosize=False,
800 width=1000,
801 height=1000,
802 margin={'l': 65, 'r': 50, 'b': 65, 't': 90})
803 fig = go.Figure(data=data, layout=layout)
805 return iplot(fig, filename='SCPS')
807 else:
808 # use matplotlib
809 scps = self._calc_shifted_cross_power_spectrum()
810 PLT.subplot_3dsurface(scps.astype(np.float32))
812 def _get_opt_winpos_winsize(self) -> None:
813 """Calculate optimal window position and size in reference image units.
815 NOTE: The returned values are computed according to DGM, cloud_mask and trueCornerLonLat.
816 """
817 # dummy algorithm: get center position of overlap instead of searching ideal window position in whole overlap
818 # TODO automatischer Algorithmus zur Bestimmung der optimalen Window Position
820 wp = tuple(self.win_pos_XY)
821 assert type(self.win_pos_XY) in [tuple, list, np.ndarray], \
822 'The window position must be a tuple of two elements. Got %s with %s elements.' % (type(wp), len(wp))
823 wp = tuple(wp)
825 if None in wp:
826 # use centroid point if possible
827 overlap_center_pos_x, overlap_center_pos_y = self.overlap_poly.centroid.coords.xy
828 wp = (wp[0] if wp[0] else overlap_center_pos_x[0]), (wp[1] if wp[1] else overlap_center_pos_y[0])
830 # validate window position
831 if not self.overlap_poly.buffer(1e-5).contains(Point(wp)):
832 # in case the centroid point is not within overlap area
833 if not self.q:
834 warnings.warn("The centroid point of the two input images could not be used as matching window "
835 "position since it is outside of the overlap area. Instead the so called "
836 "'representative point' is used. Alternatively you can provide your own window "
837 "position as input parameter.")
839 # -> use representative point: a point that is garanteed to be within overlap polygon
840 overlap_center_pos_x, overlap_center_pos_y = self.overlap_poly.representative_point().coords.xy
841 wp = overlap_center_pos_x[0], overlap_center_pos_y[0]
843 assert self.overlap_poly.buffer(1e-5).contains(Point(wp))
845 else:
846 # validate window position
847 if not self.overlap_poly.buffer(1e-5).contains(Point(wp)):
848 self._handle_error(ValueError('The provided window position %s/%s is outside of the overlap '
849 'area of the two input images. Check the coordinates.' % wp))
851 # check if window position is within bad data area if a respective mask has been provided
852 for im in [self.ref, self.shift]:
853 if im.mask_baddata is not None:
854 imX, imY = mapXY2imXY(wp, im.mask_baddata.gt)
856 if im.mask_baddata[int(imY), int(imX)] is True:
857 self._handle_error(
858 RuntimeError('According to the provided bad data mask for the %s the chosen window position '
859 '%s / %s is within a bad data area. Using this window position for coregistration '
860 'is not reasonable. Please provide a better window position!'
861 % (im.imName, wp[0], wp[1])))
863 self.win_pos_XY = wp
864 self.win_size_XY = (int(self.win_size_XY[0]), int(self.win_size_XY[1])) if self.win_size_XY else (512, 512)
866 def _get_clip_window_properties(self) -> None:
867 """Calculate all properties of the matching window and the other window.
869 These windows are used to read the corresponding image positions in the reference and the target image.
871 NOTE: Even if X- and Y-dimension of the target window is equal, the output window can be NON-quadratic!
872 """
873 # FIXME image sizes like 10000*256 are still possible
875 wpX, wpY = self.win_pos_XY
876 wsX, wsY = self.win_size_XY
878 # image units -> map units
879 ref_wsX = wsX * self.ref.xgsd
880 ref_wsY = wsY * self.ref.ygsd
881 shift_wsX = wsX * self.shift.xgsd
882 shift_wsY = wsY * self.shift.ygsd
884 ref_box_kwargs = \
885 dict(wp=(wpX, wpY),
886 ws=(ref_wsX, ref_wsY),
887 gt=self.ref.gt)
888 shift_box_kwargs = \
889 dict(wp=(wpX, wpY),
890 ws=(shift_wsX, shift_wsY),
891 gt=self.shift.gt)
892 matchBox =\
893 boxObj(**ref_box_kwargs) if self.grid2use == 'ref' else \
894 boxObj(**shift_box_kwargs)
895 otherBox = \
896 boxObj(**shift_box_kwargs) if self.grid2use == 'ref' else \
897 boxObj(**ref_box_kwargs)
898 overlapWin = \
899 boxObj(mapPoly=self.overlap_poly,
900 gt=self.ref.gt)
902 # clip matching window to overlap area
903 matchBox.mapPoly = matchBox.mapPoly.intersection(overlapWin.mapPoly)
905 # check if matchBox extent touches no data area of the image -> if yes: shrink it
906 overlapPoly_within_matchWin = matchBox.mapPoly.intersection(self.overlap_poly)
907 if overlapPoly_within_matchWin.area < matchBox.mapPoly.area:
908 wsX_start, wsY_start = \
909 1 if wsX >= wsY else \
910 wsX / wsY, 1 if wsY >= wsX else \
911 wsY / wsX
912 box = boxObj(**dict(wp=(wpX, wpY),
913 ws=(wsX_start, wsY_start),
914 gt=matchBox.gt))
915 while True:
916 box.buffer_imXY(1, 1)
917 if not box.mapPoly.within(overlapPoly_within_matchWin):
918 box.buffer_imXY(-1, -1)
919 matchBox = box
920 break
922 # move matching window to imref grid or im2shift grid
923 mW_rows, mW_cols = \
924 (self.ref.rows, self.ref.cols) if self.grid2use == 'ref' else \
925 (self.shift.rows, self.shift.cols)
926 matchBox.mapPoly = move_shapelyPoly_to_image_grid(matchBox.mapPoly,
927 matchBox.gt, mW_rows,
928 mW_cols,
929 'NW')
931 # check, if matchBox was moved outside of overlap_poly when moving it to the image grid
932 if not matchBox.mapPoly.within(overlapWin.mapPoly):
933 # further shrink matchPoly (1 px buffer is enough because the window was only moved to the grid)
934 xLarger, yLarger = matchBox.is_larger_DimXY(overlapWin.boundsIm)
935 matchBox.buffer_imXY(-1 if xLarger else 0,
936 -1 if yLarger else 0)
938 # matching_win directly on grid2use (fix rounding error through coordinate transformation)
939 matchBox.imPoly = round_shapelyPoly_coords(matchBox.imPoly, precision=0)
941 # check if matching window larger than the other one or equal
942 if not (matchBox.mapPoly.within(otherBox.mapPoly) or
943 matchBox.mapPoly == otherBox.mapPoly):
944 # if yes, find the smallest 'other window' that encloses the matching window
945 otherBox.boxImYX = \
946 get_smallest_boxImYX_that_contains_boxMapYX(
947 matchBox.boxMapYX,
948 otherBox.gt,
949 tolerance_ndigits=5 # avoids float coordinate rounding issues
950 )
952 # in case after enlarging the 'other window', it gets too large for the overlap area
953 # -> shrink match window and recompute smallest possible other window until everything is fine
954 t_start = time.time()
955 while not otherBox.mapPoly.within(overlapWin.mapPoly):
956 xLarger, yLarger = otherBox.is_larger_DimXY(overlapWin.boundsIm)
957 matchBox.buffer_imXY(-1 if xLarger else 0,
958 -1 if yLarger else 0)
959 previous_area = otherBox.mapPoly.area
960 otherBox.boxImYX = \
961 get_smallest_boxImYX_that_contains_boxMapYX(
962 matchBox.boxMapYX,
963 otherBox.gt,
964 tolerance_ndigits=5 # avoids float coordinate rounding issues
965 )
967 if previous_area == otherBox.mapPoly.area or \
968 time.time() - t_start > 1.5:
969 # happens, e.g, in case of a triangular footprint
970 # NOTE: first condition is not always fulfilled -> therefore added timeout of 1.5 sec
971 self._handle_error(
972 RuntimeError('Matching window in target image is larger than overlap area but further shrinking '
973 'the matching window is not possible. Check if the footprints of the input data have '
974 'been computed correctly.' +
975 (' Matching window shrinking timed out.' if time.time() - t_start > 5 else '')))
976 break # break out of while loop in order to avoid that code gets stuck here
978 # output validation
979 for winBox in [matchBox, otherBox]:
980 if winBox.imDimsYX[0] < 16 or \
981 winBox.imDimsYX[1] < 16:
982 self._handle_error(
983 RuntimeError("One of the input images does not have sufficient gray value information "
984 "(non-no-data values) for placing a matching window at the position %s. "
985 "Matching failed." % str((wpX, wpY))))
987 if self.success is not False:
988 # check result -> ProgrammingError if not fulfilled
989 def within_equal(inner, outer):
990 return inner.within(outer.buffer(1e-5)) or \
991 inner.equals(outer)
993 assert within_equal(matchBox.mapPoly,
994 otherBox.mapPoly)
995 assert within_equal(otherBox.mapPoly,
996 overlapWin.mapPoly)
998 if self.grid2use == 'ref':
999 self.imfft_xgsd = self.ref.xgsd
1000 self.imfft_ygsd = self.ref.ygsd
1001 self.ref.win = matchBox
1002 self.shift.win = otherBox
1003 else:
1004 self.imfft_xgsd = self.shift.xgsd
1005 self.imfft_ygsd = self.shift.ygsd
1006 self.ref.win = otherBox
1007 self.shift.win = matchBox
1009 self.matchBox = matchBox
1010 self.otherBox = otherBox
1012 self.ref.win.size_YX = tuple([int(i) for i in self.ref.win.imDimsYX])
1013 self.shift.win.size_YX = tuple([int(i) for i in self.shift.win.imDimsYX])
1014 match_win_size_XY = tuple(reversed([int(i) for i in matchBox.imDimsYX]))
1016 if not self.q and \
1017 match_win_size_XY != self.win_size_XY:
1018 print('Target window size %s not possible due to too small overlap area or window position too close '
1019 'to an image edge. New matching window size: %s.' % (self.win_size_XY, match_win_size_XY))
1021 # write_shp('matchMapPoly.shp', matchBox.mapPoly,matchBox.prj)
1022 # write_shp('otherMapPoly.shp', otherBox.mapPoly,otherBox.prj)
1024 def _get_image_windows_to_match(self) -> None:
1025 """Read the matching window and the other window as subsets.
1027 The other window is resampled to the resolution and the pixel grid of the matching window.
1028 The result consists of two images with the same dimensions and exactly the same corner coordinates.
1029 """
1030 match_fullGeoArr = self.ref if self.grid2use == 'ref' else self.shift
1031 other_fullGeoArr = self.shift if self.grid2use == 'ref' else self.ref
1033 # read matchWin via subset-read -> self.matchWin.data
1034 rS, rE, cS, cE = GEO.get_GeoArrayPosition_from_boxImYX(self.matchBox.boxImYX)
1035 assert np.array_equal(np.abs(np.array([rS, rE, cS, cE])), np.array([rS, rE, cS, cE])) and \
1036 rE <= match_fullGeoArr.rows and cE <= match_fullGeoArr.cols, \
1037 'Requested area is not completely within the input array for %s.' % match_fullGeoArr.imName
1038 self.matchWin = GeoArray(match_fullGeoArr[rS:rE + 1, cS:cE + 1, match_fullGeoArr.band4match],
1039 geotransform=GEO.get_subset_GeoTransform(match_fullGeoArr.gt, self.matchBox.boxImYX),
1040 projection=copy(match_fullGeoArr.prj),
1041 nodata=copy(match_fullGeoArr.nodata))
1042 self.matchWin.imID = match_fullGeoArr.imID
1044 # read otherWin via subset-read
1045 rS, rE, cS, cE = GEO.get_GeoArrayPosition_from_boxImYX(self.otherBox.boxImYX)
1046 assert np.array_equal(np.abs(np.array([rS, rE, cS, cE])), np.array([rS, rE, cS, cE])) and \
1047 rE <= other_fullGeoArr.rows and cE <= other_fullGeoArr.cols, \
1048 'Requested area is not completely within the input array for %s.' % other_fullGeoArr.imName
1049 self.otherWin = GeoArray(other_fullGeoArr[rS:rE + 1, cS:cE + 1, other_fullGeoArr.band4match],
1050 geotransform=GEO.get_subset_GeoTransform(other_fullGeoArr.gt, self.otherBox.boxImYX),
1051 projection=copy(other_fullGeoArr.prj),
1052 nodata=copy(other_fullGeoArr.nodata))
1053 self.otherWin.imID = other_fullGeoArr.imID
1055 # self.matchWin.deepcopy_array()
1056 # self.otherWin.deepcopy_array()
1058 if self.v:
1059 print('Original matching windows:')
1060 ref_data, shift_data = (self.matchWin[:], self.otherWin[:]) if self.grid2use == 'ref' else \
1061 (self.otherWin[:], self.matchWin[:])
1062 PLT.subplot_imshow([ref_data, shift_data], [self.ref.title, self.shift.title], grid=True)
1064 # resample otherWin.arr to the resolution of matchWin AND make sure the pixel edges are identical
1065 # (in order to make each image show the same window with the same coordinates)
1066 # TODO replace cubic resampling by PSF resampling - average resampling leads to sinus like distortions in the
1067 # TODO fft image that make a precise coregistration impossible. Thats why there is currently no way around
1068 # TODO cubic resampling.
1069 tgt_xmin, tgt_xmax, tgt_ymin, tgt_ymax = self.matchBox.boundsMap
1071 # equalize pixel grids and projection of matchWin and otherWin (ONLY if grids are really different)
1072 if not (self.matchWin.xygrid_specs == self.otherWin.xygrid_specs and
1073 prj_equal(self.matchWin.prj, self.otherWin.prj)):
1074 self.otherWin.arr, self.otherWin.gt = warp_ndarray(self.otherWin.arr,
1075 self.otherWin.gt,
1076 self.otherWin.prj,
1077 self.matchWin.prj,
1078 out_gsd=(self.imfft_xgsd, abs(self.imfft_ygsd)),
1079 out_bounds=([tgt_xmin, tgt_ymin, tgt_xmax, tgt_ymax]),
1080 rspAlg=_dict_rspAlg_rsp_Int[self.rspAlg_calc],
1081 in_nodata=self.otherWin.nodata,
1082 CPUs=self.CPUs,
1083 progress=False)[:2]
1085 if self.matchWin.shape != self.otherWin.shape:
1086 self._handle_error(
1087 RuntimeError('Caught a possible ProgrammingError at window position %s: Bad output of '
1088 'get_image_windows_to_match. Reference image shape is %s whereas shift '
1089 'image shape is %s.' % (str(self.matchBox.wp), self.matchWin.shape, self.otherWin.shape)),
1090 warn=True)
1092 # check of odd dimensions of output images
1093 rows, cols = [i if i % 2 == 0 else i - 1 for i in self.matchWin.shape]
1094 self.matchWin.arr, self.otherWin.arr = self.matchWin.arr[:rows, :cols], self.otherWin.arr[:rows, :cols]
1095 if self.matchWin.box.imDimsYX != self.matchBox.imDimsYX:
1096 self.matchBox = self.matchWin.box # update matchBox
1097 self.otherBox = self.otherWin.box # update otherBox
1099 assert self.matchWin.arr is not None and self.otherWin.arr is not None, 'Creation of matching windows failed.'
1101 @staticmethod
1102 def _shrink_winsize_to_binarySize(win_shape_YX: tuple,
1103 target_size: tuple = None
1104 ) -> Optional[Tuple[int, int]]:
1105 """Shrink a given window size to the closest binary window size (a power of 2).
1107 NOTE: X- and Y-dimension are handled separately.
1109 :param win_shape_YX: source window shape as pixel units (rows,colums)
1110 :param target_size: source window shape as pixel units (rows,colums)
1111 """
1112 binarySizes = [2 ** i for i in range(3, 14)] # [8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192]
1113 possibSizes_X = [i for i in binarySizes if i <= win_shape_YX[1]]
1114 possibSizes_Y = [i for i in binarySizes if i <= win_shape_YX[0]]
1115 if possibSizes_X and possibSizes_Y:
1116 tgt_size_X, tgt_size_Y = target_size if target_size else (max(possibSizes_X), max(possibSizes_Y))
1117 closest_to_target_X = int(min(possibSizes_X, key=lambda x: abs(x - tgt_size_X)))
1118 closest_to_target_Y = int(min(possibSizes_Y, key=lambda y: abs(y - tgt_size_Y)))
1119 return closest_to_target_Y, closest_to_target_X
1120 else:
1121 return None
1123 def _calc_shifted_cross_power_spectrum(self,
1124 im0=None,
1125 im1=None,
1126 precision=np.complex64
1127 ) -> np.ndarray:
1128 """Calculate the shifted cross power spectrum for quantifying X/Y-shifts.
1130 :param im0: reference image
1131 :param im1: subject image to shift
1132 :param precision: to be quantified as a datatype
1133 :return: 2D-numpy-array of the shifted cross power spectrum
1134 """
1135 im0 = im0 if im0 is not None else self.matchWin[:] if self.matchWin.imID == 'ref' else self.otherWin[:]
1136 im1 = im1 if im1 is not None else self.otherWin[:] if self.otherWin.imID == 'shift' else self.matchWin[:]
1138 assert im0.shape == im1.shape, 'The reference and the target image must have the same dimensions.'
1139 if im0.shape[0] % 2 != 0:
1140 warnings.warn('Odd row count in one of the match images!')
1141 if im1.shape[1] % 2 != 0:
1142 warnings.warn('Odd column count in one of the match images!')
1144 wsYX = self._shrink_winsize_to_binarySize(im0.shape) if self.bin_ws else im0.shape
1145 wsYX = ((min(wsYX),) * 2 if self.force_quadratic_win else wsYX) if wsYX else None
1147 if wsYX not in [None, (0, 0)]:
1148 time0 = time.time()
1149 if self.v:
1150 print('final window size: %s/%s (X/Y)' % (wsYX[1], wsYX[0]))
1151 # FIXME size of self.matchWin is not updated
1152 # FIXME CoRegPoints_grid.WIN_SZ is taken from self.matchBox.imDimsYX but this is not updated
1154 center_YX = np.array(im0.shape) / 2
1155 xmin, xmax = int(center_YX[1] - wsYX[1] / 2), int(center_YX[1] + wsYX[1] / 2)
1156 ymin, ymax = int(center_YX[0] - wsYX[0] / 2), int(center_YX[0] + wsYX[0] / 2)
1158 in_arr0 = im0[ymin:ymax, xmin:xmax].astype(precision)
1159 in_arr1 = im1[ymin:ymax, xmin:xmax].astype(precision)
1161 if self.v:
1162 PLT.subplot_imshow([np.real(in_arr0).astype(np.float32), np.real(in_arr1).astype(np.float32)],
1163 ['FFTin ' + self.ref.title, 'FFTin ' + self.shift.title], grid=True)
1165 fft_arr0, fft_arr1 = None, None
1166 if pyfftw and self.fftw_works is not False: # if module is installed and working
1167 try:
1168 fft_arr0 = pyfftw.FFTW(in_arr0, np.empty_like(in_arr0), axes=(0, 1))()
1169 fft_arr1 = pyfftw.FFTW(in_arr1, np.empty_like(in_arr1), axes=(0, 1))()
1171 # catch empty output arrays (for some reason this happens sometimes) -> use numpy fft
1172 # => this is caused by the call of pyfftw.FFTW. Exactly at that moment the input array
1173 # in_arr0 is overwritten with zeros (maybe this is a bug in pyFFTW?)
1174 if np.std(fft_arr0) == 0 or np.std(fft_arr1) == 0:
1175 raise RuntimeError('FFTW result is unexpectedly empty.')
1177 self.fftw_works = True
1179 except RuntimeError:
1180 self.fftw_works = False
1182 # recreate input arrays and use numpy fft as fallback
1183 in_arr0 = im0[ymin:ymax, xmin:xmax].astype(precision)
1184 in_arr1 = im1[ymin:ymax, xmin:xmax].astype(precision)
1186 if self.fftw_works is False or fft_arr0 is None or fft_arr1 is None:
1187 fft_arr0 = fft2(in_arr0)
1188 fft_arr1 = fft2(in_arr1)
1190 # GeoArray(fft_arr0.astype(np.float32)).show(figsize=(15,15))
1191 # GeoArray(fft_arr1.astype(np.float32)).show(figsize=(15,15))
1193 if self.v:
1194 print('forward FFTW: %.2fs' % (time.time() - time0))
1196 eps = np.abs(fft_arr1).max() * 1e-15
1197 # cps == cross-power spectrum of im0 and im2
1199 with np.errstate(invalid='ignore'): # ignore RuntimeWarning: invalid value encountered in true_divide
1200 temp = np.array(fft_arr0 * fft_arr1.conjugate()) / (np.abs(fft_arr0) * np.abs(fft_arr1) + eps)
1202 time0 = time.time()
1203 if 'pyfft' in globals():
1204 ifft_arr = pyfftw.FFTW(temp, np.empty_like(temp), axes=(0, 1), direction='FFTW_BACKWARD')()
1205 else:
1206 ifft_arr = ifft2(temp)
1207 if self.v:
1208 print('backward FFTW: %.2fs' % (time.time() - time0))
1210 cps = np.abs(ifft_arr)
1211 # scps = shifted cps => shift the zero-frequency component to the center of the spectrum
1212 scps = fftshift(cps)
1213 if self.v:
1214 PLT.subplot_imshow([np.real(in_arr0).astype(np.uint16), np.real(in_arr1).astype(np.uint16),
1215 np.real(fft_arr0).astype(np.uint8), np.real(fft_arr1).astype(np.uint8), scps],
1216 titles=['matching window im0', 'matching window im1',
1217 "fft result im0", "fft result im1", "cross power spectrum"], grid=True)
1218 PLT.subplot_3dsurface(np.real(scps).astype(np.float32))
1219 else:
1220 scps = None
1221 self._handle_error(
1222 RuntimeError('The matching window became too small for calculating a reliable match. Matching failed.'))
1224 self.fftw_win_size_YX = wsYX
1225 return scps
1227 @staticmethod
1228 def _get_peakpos(scps: np.ndarray) -> np.ndarray:
1229 """Return the row/column position of the peak within the given cross power spectrum.
1231 :param scps: shifted cross power spectrum
1232 :return: [row, column]
1233 """
1234 max_flat_idx = np.argmax(scps)
1235 return np.array(np.unravel_index(max_flat_idx, scps.shape))
1237 @staticmethod
1238 def _get_shifts_from_peakpos(peakpos: np.array,
1239 arr_shape: tuple
1240 ) -> (float, float):
1241 y_shift = peakpos[0] - arr_shape[0] // 2
1242 x_shift = peakpos[1] - arr_shape[1] // 2
1243 return x_shift, y_shift
1245 @staticmethod
1246 def _clip_image(im, center_YX, winSzYX): # TODO this is also implemented in GeoArray
1248 def get_bounds(YX: tuple, wsY: float, wsX: float):
1249 return (
1250 int(YX[1] - (wsX / 2)),
1251 int(YX[1] + (wsX / 2)),
1252 int(YX[0] - (wsY / 2)),
1253 int(YX[0] + (wsY / 2))
1254 )
1256 wsY, wsX = winSzYX
1257 xmin, xmax, ymin, ymax = get_bounds(center_YX, wsY, wsX)
1259 return im[ymin:ymax, xmin:xmax]
1261 def _get_grossly_deshifted_images(self, im0, im1, x_intshift, y_intshift):
1262 # TODO this is also implemented in GeoArray # this should update ref.win.data and shift.win.data
1263 # FIXME avoid that matching window gets smaller although shifting it with the previous win_size would not move
1264 # it into nodata-area
1265 # get_grossly_deshifted_im0
1266 old_center_YX = np.array(im0.shape) / 2
1267 new_center_YX = [old_center_YX[0] + y_intshift,
1268 old_center_YX[1] + x_intshift]
1270 x_left = new_center_YX[1]
1271 x_right = im0.shape[1] - new_center_YX[1]
1272 y_above = new_center_YX[0]
1273 y_below = im0.shape[0] - new_center_YX[0]
1274 maxposs_winsz_x = 2 * min(x_left, x_right)
1275 maxposs_winsz_y = 2 * min(y_above, y_below)
1277 if self.force_quadratic_win:
1278 maxposs_winsz_x = maxposs_winsz_y = min([maxposs_winsz_x, maxposs_winsz_y])
1280 gdsh_im0 = self._clip_image(im0, new_center_YX, [maxposs_winsz_y, maxposs_winsz_x])
1282 # get_corresponding_im1_clip
1283 crsp_im1 = self._clip_image(im1, np.array(im1.shape) / 2, gdsh_im0.shape)
1285 if self.v:
1286 PLT.subplot_imshow([self._clip_image(im0, old_center_YX, gdsh_im0.shape), crsp_im1],
1287 titles=['reference original', 'target'],
1288 grid=True)
1289 PLT.subplot_imshow([gdsh_im0, crsp_im1],
1290 titles=['reference virtually shifted', 'target'],
1291 grid=True)
1293 return gdsh_im0, crsp_im1
1295 def _find_side_maximum(self, scps):
1296 # peakR, peakC = self._get_peakpos(scps)
1297 peakR, peakC = scps.shape[0] // 2, scps.shape[1] // 2
1298 profileX = scps[peakR, :].flatten() # row profile with values from left to right
1299 profileY = scps[:, peakC].flatten() # column profile with values from top to bottom
1301 # get scps values of side maxima
1302 sm_left, sm_right = profileX[peakC - 1], profileX[peakC + 1]
1303 sm_above, sm_below = profileY[peakR - 1], profileY[peakR + 1]
1305 sidemax_lr = {'value': max([sm_left, sm_right]),
1306 'side': 'left' if sm_left > sm_right else 'right',
1307 'direction_factor': -1 if sm_left > sm_right else 1}
1308 sidemax_ab = {'value': max([sm_above, sm_below]),
1309 'side': 'above' if sm_above > sm_below else 'below',
1310 'direction_factor': -1 if sm_above > sm_below else 1}
1312 if self.v:
1313 print('Horizontal side maximum found %s. value: %s' % (sidemax_lr['side'], sidemax_lr['value']))
1314 print('Vertical side maximum found %s. value: %s' % (sidemax_ab['side'], sidemax_ab['value']))
1315 PLT.subplot_2dline([[range(profileX.size), profileX], [range(profileY.size), profileY]],
1316 titles=['X-Profile', 'Y-Profile'], shapetuple=(1, 2), grid=True)
1318 return sidemax_lr, sidemax_ab
1320 def _calc_integer_shifts(self, scps: np.ndarray) -> (int, int):
1321 peakpos = self._get_peakpos(scps)
1322 x_intshift, y_intshift = self._get_shifts_from_peakpos(peakpos, scps.shape)
1324 return x_intshift, y_intshift
1326 def _calc_shift_reliability(self, scps: np.ndarray):
1327 """Calculate a confidence percentage to be used as an assessment for reliability of the calculated shifts.
1329 :param scps: shifted cross power spectrum
1330 :return:
1331 """
1332 # calculate mean power at peak
1333 peakR, peakC = self._get_peakpos(scps)
1334 power_at_peak = np.mean(scps[peakR - 1:peakR + 2, peakC - 1:peakC + 2])
1336 # calculate mean power without peak + 3* standard deviation
1337 scps_masked = scps
1338 scps_masked[peakR - 1:peakR + 2, peakC - 1:peakC + 2] = -9999
1339 scps_masked = np.ma.masked_equal(scps_masked, -9999)
1340 power_without_peak = np.mean(scps_masked) + 2 * np.std(scps_masked)
1342 # calculate confidence
1343 confid = 100 - ((power_without_peak / power_at_peak) * 100)
1344 confid = 100 if confid > 100 else 0 if confid < 0 else confid
1346 if not self.q:
1347 print('Estimated reliability of the calculated shifts: %.1f' % confid, '%')
1349 return confid
1351 def _validate_integer_shifts(self, im0, im1, x_intshift, y_intshift):
1353 if (x_intshift, y_intshift) != (0, 0):
1354 # temporalily deshift images on the basis of calculated integer shifts
1355 gdsh_im0, crsp_im1 = self._get_grossly_deshifted_images(im0, im1, x_intshift, y_intshift)
1357 # check if integer shifts are now gone (0/0)
1358 scps = self._calc_shifted_cross_power_spectrum(gdsh_im0, crsp_im1)
1360 if scps is not None:
1361 if scps.shape[0] < 3 or scps.shape[1] < 3:
1362 self._handle_error(RuntimeError('Shifted cross power spectrum became too small for computing the '
1363 'point of registration. Matching failed.'))
1364 return 'invalid', None, None, scps
1366 peakpos = self._get_peakpos(scps)
1367 x_shift, y_shift = self._get_shifts_from_peakpos(peakpos, scps.shape)
1368 if (x_shift, y_shift) == (0, 0):
1369 return 'valid', 0, 0, scps
1370 else:
1371 return 'invalid', x_shift, y_shift, scps
1372 else:
1373 return 'invalid', None, None, scps
1374 else:
1375 return 'valid', 0, 0, None
1377 def _calc_subpixel_shifts(self, scps: np.ndarray):
1378 sidemax_lr, sidemax_ab = self._find_side_maximum(scps)
1379 x_subshift = (sidemax_lr['direction_factor'] * sidemax_lr['value']) / (np.max(scps) + sidemax_lr['value'])
1380 y_subshift = (sidemax_ab['direction_factor'] * sidemax_ab['value']) / (np.max(scps) + sidemax_ab['value'])
1382 return x_subshift, y_subshift
1384 @staticmethod
1385 def _get_total_shifts(x_intshift: int,
1386 y_intshift: int,
1387 x_subshift: float,
1388 y_subshift: float):
1389 return x_intshift + x_subshift, y_intshift + y_subshift
1391 def _get_deshifted_otherWin(self) -> GeoArray:
1392 """Return a de-shifted version of self.otherWin as a GeoArray instance.
1394 The output dimensions and geographic bounds are equal to those of self.matchWin and geometric shifts are
1395 corrected according to the previously computed X/Y shifts within the matching window. This allows direct
1396 application of algorithms e.g. measuring image similarity.
1398 The image subset that is resampled in this function is always the same that has been resampled during
1399 computation of geometric shifts (usually the image with the higher geometric resolution).
1401 :returns: GeoArray instance of de-shifted self.otherWin
1402 """
1403 # shift vectors have been calculated to fit target image onto reference image
1404 # -> so the shift vectors have to be inverted if shifts are applied to reference image
1405 coreg_info = self._get_inverted_coreg_info() if self.otherWin.imID == 'ref' else self.coreg_info
1407 matchFull = self.ref if self.matchWin.imID == 'ref' else self.shift
1408 otherFull = self.ref if self.otherWin.imID == 'ref' else self.shift
1410 ds_results = DESHIFTER(otherFull, coreg_info,
1411 band2process=otherFull.band4match + 1,
1412 clipextent=self.matchBox.mapPoly.bounds,
1413 target_xyGrid=matchFull.xygrid_specs,
1414 q=True
1415 ).correct_shifts()
1416 return ds_results['GeoArray_shifted']
1418 def _validate_ssim_improvement(self,
1419 v: bool = False
1420 ) -> (float, float):
1421 """Compute mean structural similarity index between reference and target image before and after co-registration.
1423 :param v: verbose mode: shows images of the matchWin, otherWin and shifted version of otherWin
1424 :return: SSIM before an after shift correction
1425 """
1426 from skimage.metrics import structural_similarity as ssim # import here to avoid static TLS import error
1428 assert self.success is not None, \
1429 'Calculate geometric shifts first before trying to measure image similarity improvement!'
1430 assert self.success in [True, None], \
1431 'Since calculation of geometric shifts failed, no image similarity improvement can be measured.'
1433 def normalize(array: np.ndarray) -> np.ndarray:
1434 minval = np.min(array)
1435 maxval = np.max(array)
1437 # avoid zerodivision
1438 if maxval == minval:
1439 maxval += 1e-5
1441 return ((array - minval) / (maxval - minval)).astype(np.float64)
1443 # compute SSIM BEFORE shift correction #
1444 ########################################
1446 # using gaussian weights could lead to value errors in case of small images when the automatically calculated
1447 # window size exceeds the image size
1448 self.ssim_orig = ssim(normalize(np.ma.masked_equal(self.matchWin[:],
1449 self.matchWin.nodata)),
1450 normalize(np.ma.masked_equal(self.otherWin[:],
1451 self.otherWin.nodata)),
1452 data_range=1)
1454 # compute SSIM AFTER shift correction #
1455 #######################################
1457 # resample otherWin while correcting detected shifts and match geographic bounds of matchWin
1458 otherWin_deshift_geoArr = self._get_deshifted_otherWin()
1460 # get the corresponding matchWin data
1461 matchWinData = self.matchWin[:]
1463 # check if shapes of two images are unequal (due to bug (?), in some cases otherWin_deshift_geoArr does not have
1464 # the exact same dimensions as self.matchWin -> maybe bounds are handled differently by gdal.Warp)
1465 if not self.matchWin.shape == otherWin_deshift_geoArr.shape: # FIXME this seems to be already fixed
1466 warnings.warn('SSIM input array shapes are not equal. This should not happen. '
1467 'If it happens with your data, please report it here so that the issue can be fixed: '
1468 'https://git.gfz-potsdam.de/danschef/arosics/-/issues/85')
1469 matchFull = \
1470 self.ref if self.matchWin.imID == 'ref' else\
1471 self.shift
1472 matchWinData, _, _ = matchFull.get_mapPos(self.matchBox.mapPoly.bounds,
1473 self.matchWin.prj,
1474 rspAlg='cubic',
1475 band2get=matchFull.band4match)
1476 self.matchWin.clip_to_poly(self.matchBox.mapPoly)
1478 # at the image edges it is possible that the size of the matchBox must be reduced
1479 # in order to make array shapes match
1480 if not matchWinData.shape == otherWin_deshift_geoArr.shape: # FIXME this seems to be already fixed
1481 self.matchBox.buffer_imXY(float(-np.ceil(abs(self.x_shift_px))),
1482 float(-np.ceil(abs(self.y_shift_px))))
1483 otherWin_deshift_geoArr = self._get_deshifted_otherWin()
1484 matchWinData, _, _ = matchFull.get_mapPos(self.matchBox.mapPoly.bounds,
1485 self.matchWin.prj,
1486 rspAlg='cubic',
1487 band2get=matchFull.band4match)
1489 # output validation
1490 if not matchWinData.shape == otherWin_deshift_geoArr.shape:
1491 warnings.warn('SSIM input array shapes could not be equalized. SSIM calculation failed. '
1492 'SSIM of the de-shifted target image is set to 0.')
1493 self.ssim_deshifted = 0
1495 return self.ssim_orig, self.ssim_deshifted
1497 self.ssim_deshifted = ssim(normalize(np.ma.masked_equal(otherWin_deshift_geoArr[:],
1498 otherWin_deshift_geoArr.nodata)),
1499 normalize(np.ma.masked_equal(matchWinData,
1500 self.matchWin.nodata)),
1501 data_range=1)
1503 if v:
1504 GeoArray(matchWinData).show()
1505 self.otherWin.show()
1506 otherWin_deshift_geoArr.show()
1508 if not self.q:
1509 print('Image similarity within the matching window (SSIM before/after correction): %.4f => %.4f'
1510 % (self.ssim_orig, self.ssim_deshifted))
1512 self.ssim_improved = self.ssim_orig <= self.ssim_deshifted
1514 # write win data to disk
1515 # outDir = '/home/gfz-fe/scheffler/temp/ssim_debugging/'
1516 # GeoArray(matchWinData, matchWinGt, matchWinPrj).save(outDir+'matchWinData.bsq')
1518 # otherWinGt = (self.otherWin.boundsMap[0], self.matchWin.xgsd, 0,
1519 # self.otherWin.boundsMap[3], 0, -self.matchWin.imParams.ygsd)
1520 # GeoArray(self.otherWin.data, therWinGt, self.otherWin.prj).save(outDir+'otherWin.data.bsq')
1522 # otherWin_deshift_geoArr.save(outDir+''shifted.bsq')
1524 return self.ssim_orig, self.ssim_deshifted
1526 @property
1527 def ssim_improved(self) -> bool:
1528 """Return True if image similarity within the matching window has been improved by co-registration."""
1529 if self.success is True:
1530 if self._ssim_improved is None:
1531 ssim_orig, ssim_deshifted = self._validate_ssim_improvement()
1532 self._ssim_improved = ssim_orig <= ssim_deshifted
1534 return self._ssim_improved
1536 @ssim_improved.setter
1537 def ssim_improved(self, has_improved: bool):
1538 self._ssim_improved = has_improved
1540 def calculate_spatial_shifts(self) -> str:
1541 """Compute the global X/Y shift between reference and the target image within the matching window.
1543 :return: 'success' or 'fail'
1544 """
1545 if self.success is False:
1546 return 'fail'
1548 # set self.matchWin and self.otherWin (GeoArray instances)
1549 self._get_image_windows_to_match() # 45-90ms
1551 im0 = self.matchWin[:] if self.matchWin.imID == 'ref' else self.otherWin[:]
1552 im1 = self.otherWin[:] if self.otherWin.imID == 'shift' else self.matchWin[:]
1554 xgsd_factor = self.imfft_xgsd / self.shift.xgsd
1555 ygsd_factor = self.imfft_ygsd / self.shift.ygsd
1557 if self.v:
1558 print('Matching windows with equalized spatial resolution:')
1559 PLT.subplot_imshow([im0, im1], [self.ref.title, self.shift.title], grid=True)
1560 print('xgsd_factor', xgsd_factor)
1561 print('ygsd_factor', ygsd_factor)
1562 print('imfft_xgsd_mapvalues', self.imfft_xgsd)
1563 print('imfft_ygsd_mapvalues', self.imfft_ygsd)
1565 # calculate cross power spectrum without any de-shifting applied
1566 scps = self._calc_shifted_cross_power_spectrum() # 8-18ms
1568 if scps is None:
1569 self.success = False
1571 return 'fail'
1573 # calculate spatial shifts
1575 # calculate integer shifts
1576 count_iter = 1
1577 x_intshift, y_intshift = self._calc_integer_shifts(scps)
1579 if (x_intshift, y_intshift) == (0, 0):
1580 # in case integer shifts are zero
1581 self.success = True
1582 else:
1583 valid_invalid, x_val_shift, y_val_shift, scps = \
1584 self._validate_integer_shifts(im0, im1, x_intshift, y_intshift)
1586 while valid_invalid != 'valid':
1587 count_iter += 1
1589 if count_iter > self.max_iter:
1590 self._handle_error(RuntimeError('No match found in the given window.'))
1591 break
1593 if valid_invalid == 'invalid' and (x_val_shift, y_val_shift) == (None, None):
1594 # this happens if matching window became too small
1595 self.success = False
1596 break
1598 if not self.q:
1599 print('No clear match found yet. Jumping to iteration %s...' % count_iter)
1600 print('input shifts: ', x_val_shift, y_val_shift)
1602 valid_invalid, x_val_shift, y_val_shift, scps = \
1603 self._validate_integer_shifts(im0, im1, x_val_shift, y_val_shift)
1605 # overwrite previous integer shifts if a valid match has been found
1606 if valid_invalid == 'valid':
1607 self.success = True
1608 x_intshift, y_intshift = x_val_shift, y_val_shift
1610 # calculate sub-pixel shifts
1611 if self.success or self.success is None:
1612 # get total pixel shifts
1613 x_subshift, y_subshift = self._calc_subpixel_shifts(scps)
1614 x_totalshift, y_totalshift = self._get_total_shifts(x_intshift, y_intshift, x_subshift, y_subshift)
1616 if max([abs(x_totalshift), abs(y_totalshift)]) > self.max_shift:
1617 self._handle_error(
1618 RuntimeError("The calculated shift (X: %s px / Y: %s px) is recognized as too large to "
1619 "be valid. If you know that it is valid, just set the '-max_shift' "
1620 "parameter to an appropriate value. Otherwise try to use a different window "
1621 "size for matching via the '-ws' parameter or define the spectral bands "
1622 "to be used for matching manually ('-br' and '-bs')."
1623 % (x_totalshift, y_totalshift)))
1624 else:
1625 self.success = True
1626 self.x_shift_px, self.y_shift_px = x_totalshift * xgsd_factor, y_totalshift * ygsd_factor
1628 # get map shifts
1629 new_originX, new_originY = imXY2mapXY((self.x_shift_px, self.y_shift_px), self.shift.gt)
1630 self.x_shift_map, self.y_shift_map = new_originX - self.shift.gt[0], new_originY - self.shift.gt[3]
1632 # get length of shift vector in map units
1633 self.vec_length_map = float(np.sqrt(self.x_shift_map ** 2 + self.y_shift_map ** 2))
1635 # get angle of shift vector
1636 self.vec_angle_deg = GEO.angle_to_north((self.x_shift_px, self.y_shift_px)).tolist()[0]
1638 # print results
1639 if not self.q:
1640 print('Detected integer shifts (X/Y): %s/%s' % (x_intshift, y_intshift))
1641 print('Detected subpixel shifts (X/Y): %s/%s' % (x_subshift, y_subshift))
1642 print('Calculated total shifts in fft pixel units (X/Y): %s/%s'
1643 % (x_totalshift, y_totalshift))
1644 print('Calculated total shifts in reference pixel units (X/Y): %s/%s'
1645 % (x_totalshift, y_totalshift))
1646 print('Calculated total shifts in target pixel units (X/Y): %s/%s'
1647 % (self.x_shift_px, self.y_shift_px))
1648 print('Calculated map shifts (X,Y):\t\t\t\t %s/%s' % (self.x_shift_map, self.y_shift_map))
1649 print('Calculated absolute shift vector length in map units: %s' % self.vec_length_map)
1650 print('Calculated angle of shift vector in degrees from North: %s' % self.vec_angle_deg)
1652 if self.x_shift_px or self.y_shift_px:
1653 self._get_updated_map_info()
1655 # set self.ssim_before and ssim_after
1656 self._validate_ssim_improvement() # FIXME uses the not updated matchWin size
1657 self.shift_reliability = self._calc_shift_reliability(scps)
1659 return 'success'
1661 def _get_updated_map_info(self) -> None:
1662 original_map_info = geotransform2mapinfo(self.shift.gt, self.shift.prj)
1663 self.updated_map_info = copy(original_map_info)
1664 self.updated_map_info[3] = str(float(original_map_info[3]) + self.x_shift_map)
1665 self.updated_map_info[4] = str(float(original_map_info[4]) + self.y_shift_map)
1667 if not self.q:
1668 print('Original map info:', original_map_info)
1669 print('Updated map info: ', self.updated_map_info)
1671 @property
1672 def coreg_info(self) -> dict:
1673 """Return a dictionary containing everything to correct the detected global X/Y shift of the target image."""
1674 if self._coreg_info:
1675 return self._coreg_info
1677 else:
1678 if self.success is None:
1679 self.calculate_spatial_shifts()
1681 self._coreg_info = {
1682 'corrected_shifts_px': {
1683 'x': self.x_shift_px,
1684 'y': self.y_shift_px
1685 },
1686 'corrected_shifts_map': {
1687 'x': self.x_shift_map,
1688 'y': self.y_shift_map
1689 },
1690 'original map info':
1691 geotransform2mapinfo(self.shift.gt, self.shift.prj),
1692 'updated map info':
1693 self.updated_map_info,
1694 'reference projection':
1695 self.ref.prj,
1696 'reference geotransform':
1697 self.ref.gt,
1698 'reference grid': [
1699 [self.ref.gt[0], self.ref.gt[0] + self.ref.gt[1]],
1700 [self.ref.gt[3], self.ref.gt[3] + self.ref.gt[5]]
1701 ],
1702 'reference extent': {
1703 'cols': self.ref.xgsd,
1704 'rows': self.ref.ygsd
1705 }, # FIXME not needed anymore
1706 'success':
1707 self.success}
1709 return self._coreg_info
1711 def _get_inverted_coreg_info(self) -> dict:
1712 """Return an inverted dictionary of coreg_info.
1714 This dictionary can be passed to DESHIFTER in order to fit the REFERENCE image onto the TARGET image.
1715 """
1716 inv_coreg_info = copy(self.coreg_info)
1717 inv_coreg_info['corrected_shifts_px']['x'] *= -1
1718 inv_coreg_info['corrected_shifts_px']['y'] *= -1
1719 inv_coreg_info['corrected_shifts_map']['x'] *= -1
1720 inv_coreg_info['corrected_shifts_map']['y'] *= -1
1721 inv_coreg_info['original map info'] = geotransform2mapinfo(self.ref.gt, self.ref.prj)
1722 inv_coreg_info['reference geotransform'] = self.shift.gt
1723 inv_coreg_info['reference grid'] = self.shift.xygrid_specs
1724 inv_coreg_info['reference projection'] = self.shift.prj
1726 if inv_coreg_info['updated map info']:
1727 updated_map_info = copy(inv_coreg_info['original map info'])
1728 updated_map_info[3] = str(float(inv_coreg_info['original map info'][3]) - self.x_shift_map)
1729 updated_map_info[4] = str(float(inv_coreg_info['original map info'][4]) - self.y_shift_map)
1730 inv_coreg_info['updated map info'] = updated_map_info
1732 return inv_coreg_info
1734 def correct_shifts(self) -> dict:
1735 """Correct the already calculated X/Y shift of the target image.
1737 :return: COREG.deshift_results (dictionary)
1738 """
1739 DS = DESHIFTER(self.shift, self.coreg_info,
1740 path_out=self.path_out,
1741 fmt_out=self.fmt_out,
1742 out_crea_options=self.out_creaOpt,
1743 out_gsd=self.out_gsd,
1744 resamp_alg=self.rspAlg_DS,
1745 align_grids=self.align_grids,
1746 match_gsd=self.match_gsd,
1747 target_xyGrid=self.target_xyGrid,
1748 nodata=self.shift.nodata,
1749 CPUs=self.CPUs,
1750 progress=self.progress,
1751 v=self.v,
1752 q=self.q)
1754 self.deshift_results = DS.correct_shifts()
1756 return self.deshift_results