Coverage for arosics/CoReg_local.py: 79%

1# -*- coding: utf-8 -*- 

2 

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. 

25 

26import warnings 

27import os 

28from copy import copy 

29from typing import Tuple, Union, Optional 

30from collections import OrderedDict 

31from multiprocessing import cpu_count 

32 

33# custom 

34from osgeo import gdal # noqa 

35from packaging.version import parse as parse_version 

36try: 

37 import pyfftw 

38 # pyfftw>=0.13.0 is currently not used due to https://github.com/pyFFTW/pyFFTW/issues/294 

39 if parse_version(pyfftw.__version__) >= parse_version('0.13.0'): 

40 pyfftw = None 

41except ImportError: 

42 pyfftw = None 

43import numpy as np 

44 

45from matplotlib import pyplot as plt # noqa F401 

46from geopandas import GeoDataFrame # noqa F401 

47 

48from .Tie_Point_Grid import Tie_Point_Grid 

49from .CoReg import COREG 

50from .DeShifter import DESHIFTER 

51from .geometry import has_metaRotation, remove_metaRotation 

52from py_tools_ds.geo.coord_trafo import transform_any_prj, reproject_shapelyGeometry 

53from py_tools_ds.geo.map_info import geotransform2mapinfo 

54from geoarray import GeoArray 

55 

56__author__ = 'Daniel Scheffler' 

57 

58 

59class COREG_LOCAL(object): 

60 """ 

61 COREG_LOCAL applies the algorithm to detect spatial shifts to the whole overlap area of the input images. 

62 

63 Spatial shifts are calculated for each point in grid of which the parameters can be adjusted using keyword 

64 arguments. Shift correction performs a polynomial transformation using the calculated shifts of each point in the 

65 grid as GCPs. Thus, this class can be used to correct for locally varying geometric distortions of the target image. 

66 

67 See help(COREG_LOCAL) for documentation. 

68 """ 

69 

70 def __init__(self, 

71 im_ref: Union[GeoArray, str], 

72 im_tgt: Union[GeoArray, str], 

73 grid_res: float, 

74 max_points: int = None, 

75 window_size: Tuple[int, int] = (256, 256), 

76 path_out: str = None, 

77 fmt_out: str = 'ENVI', 

78 out_crea_options: list = None, 

79 projectDir: str = None, 

80 r_b4match: int = 1, 

81 s_b4match: int = 1, 

82 max_iter: int = 5, 

83 max_shift: int = 5, 

84 tieP_filter_level: int = 3, 

85 min_reliability: float = 60, 

86 rs_max_outlier: float = 10, 

87 rs_tolerance: float = 2.5, 

88 rs_random_state: Optional[int] = 0, 

89 align_grids: bool = True, 

90 match_gsd: bool = False, 

91 out_gsd: float = None, 

92 target_xyGrid=None, 

93 resamp_alg_deshift: str = 'cubic', 

94 resamp_alg_calc: str = 'cubic', 

95 footprint_poly_ref: str = None, 

96 footprint_poly_tgt: str = None, 

97 data_corners_ref: list = None, 

98 data_corners_tgt: list = None, 

99 outFillVal: int = -9999, 

100 nodata: Tuple[int, int] = (None, None), 

101 calc_corners: bool = True, 

102 binary_ws: bool = True, 

103 force_quadratic_win: bool = True, 

104 mask_baddata_ref: Union[GeoArray, str] = None, 

105 mask_baddata_tgt: Union[GeoArray, str] = None, 

106 CPUs: int = None, 

107 progress: bool = True, 

108 v: bool = False, 

109 q: bool = False, 

110 ignore_errors: bool = True 

111 ) -> None: 

112 """ 

113 Get an instance of COREG_LOCAL. 

114 

115 :param im_ref: 

116 source path of reference image (any GDAL compatible image format is supported) 

117 

118 :param im_tgt: 

119 source path of image to be shifted (any GDAL compatible image format is supported) 

120 

121 :param grid_res: 

122 tie point grid resolution in pixels of the target image (x-direction) 

123 

124 :param max_points: 

125 maximum number of points used to find coregistration tie points 

126 

127 NOTE: Points are selected randomly from the given point grid (specified by 'grid_res'). 

128 If the point grid does not provide enough points, all available points are chosen. 

129 

130 :param window_size: 

131 custom matching window size [pixels] (default: (256,256)) 

132 

133 :param path_out: 

134 target path of the coregistered image 

135 - if None (default), no output is written to disk 

136 - if 'auto': /dir/of/im1/<im1>__shifted_to__<im0>.bsq 

137 

138 :param fmt_out: 

139 raster file format for output file. ignored if path_out is None. Can be any GDAL compatible raster file 

140 format (e.g. 'ENVI', 'GTIFF'; default: ENVI). Refer to https://gdal.org/drivers/raster/index.html to get a 

141 full list of supported formats. 

142 

143 :param out_crea_options: 

144 GDAL creation options for the output image, e.g. ["QUALITY=80", "REVERSIBLE=YES", "WRITE_METADATA=YES"] 

145 

146 :param projectDir: 

147 name of a project directory where to store all the output results. If given, name is inserted into all 

148 automatically generated output paths. 

149 

150 :param r_b4match: 

151 band of reference image to be used for matching (starts with 1; default: 1) 

152 

153 :param s_b4match: 

154 band of shift image to be used for matching (starts with 1; default: 1) 

155 

156 :param max_iter: 

157 maximum number of iterations for matching (default: 5) 

158 

159 :param max_shift: 

160 maximum shift distance in reference image pixel units (default: 5 px) 

161 

162 :param tieP_filter_level: 

163 filter tie points used for shift correction in different levels (default: 3). 

164 NOTE: lower levels are also included if a higher level is chosen 

165 

166 - Level 0: no tie point filtering 

167 - Level 1: Reliablity filtering 

168 - filter all tie points out that have a low reliability according to internal tests 

169 - Level 2: SSIM filtering 

170 - filters all tie points out where shift correction does not increase image similarity within 

171 matching window (measured by mean structural similarity index) 

172 - Level 3: RANSAC outlier detection 

173 

174 :param min_reliability: 

175 Tie point filtering: minimum reliability threshold, below which tie points are marked as false-positives 

176 (default: 60%) 

177 - accepts values between 0% (no reliability) and 100 % (perfect reliability) 

178 HINT: decrease this value in case of poor signal-to-noise ratio of your input data 

179 

180 :param rs_max_outlier: 

181 RANSAC tie point filtering: proportion of expected outliers (default: 10%) 

182 

183 :param rs_tolerance: 

184 RANSAC tie point filtering: percentage tolerance for max_outlier_percentage (default: 2.5%) 

185 

186 :param rs_random_state: 

187 RANSAC random state (an integer corresponds to a fixed/pseudo-random state, None randomizes the result) 

188 

189 :param align_grids: 

190 True: align the input coordinate grid to the reference (does not affect the output pixel size as long as 

191 input and output pixel sizes are compatible (5:30 or 10:30 but not 4:30)), default = True 

192 

193 :param match_gsd: 

194 True: match the input pixel size to the reference pixel size, 

195 default = False 

196 

197 :param out_gsd: 

198 output pixel size in units of the reference coordinate system (default = pixel size of the input array), 

199 given values are overridden by match_gsd=True 

200 

201 :param target_xyGrid: 

202 a list with a target x-grid and a target y-grid like [[15,45], [15,45]] 

203 This overrides 'out_gsd', 'align_grids' and 'match_gsd'. 

204 

205 :param resamp_alg_deshift: 

206 the resampling algorithm to be used for shift correction (if neccessary) 

207 valid algorithms: nearest, bilinear, cubic, cubic_spline, lanczos, average, mode, max, min, med, q1, q3 

208 (default: cubic) 

209 

210 :param resamp_alg_calc: 

211 the resampling algorithm to be used for all warping processes during calculation of spatial shifts 

212 valid algorithms: nearest, bilinear, cubic, cubic_spline, lanczos, average, mode, max, min, med, q1, q3 

213 (default: cubic (highly recommended)) 

214 

215 :param footprint_poly_ref: 

216 footprint polygon of the reference image (WKT string or shapely.geometry.Polygon), 

217 e.g. 'POLYGON ((299999 6000000, 299999 5890200, 409799 5890200, 409799 6000000, 299999 6000000))' 

218 

219 :param footprint_poly_tgt: 

220 footprint polygon of the image to be shifted (WKT string or shapely.geometry.Polygon) 

221 e.g. 'POLYGON ((299999 6000000, 299999 5890200, 409799 5890200, 409799 6000000, 299999 6000000))' 

222 

223 :param data_corners_ref: 

224 map coordinates of data corners within reference image. ignored if footprint_poly_ref is given. 

225 

226 :param data_corners_tgt: 

227 map coordinates of data corners within image to be shifted. ignored if footprint_poly_tgt is given. 

228 

229 :param outFillVal: 

230 if given the generated tie point grid is filled with this value in case no match could be found during 

231 co-registration (default: -9999) 

232 

233 :param nodata: 

234 no data values for reference image and image to be shifted 

235 

236 :param calc_corners: 

237 calculate true positions of the dataset corners in order to get a useful matching window position within 

238 the actual image overlap 

239 (default: True; deactivated if 'data_corners_im0' and 'data_corners_im1' are given) 

240 

241 :param binary_ws: 

242 use binary X/Y dimensions for the matching window (default: True) 

243 

244 :param force_quadratic_win: 

245 force a quadratic matching window (default: True) 

246 

247 :param mask_baddata_ref: 

248 path to a 2D boolean mask file (or an instance of BadDataMask) for the reference image where all bad data 

249 pixels (e.g. clouds) are marked with True and the remaining pixels with False. Must have the same 

250 geographic extent and projection like 'im_ref'. The mask is used to check if the chosen matching window 

251 position is valid in the sense of useful data. Otherwise, this window position is rejected. 

252 

253 :param mask_baddata_tgt: 

254 path to a 2D boolean mask file (or an instance of BadDataMask) for the image to be shifted where all bad 

255 data pixels (e.g. clouds) are marked with True and the remaining pixels with False. Must have the same 

256 geographic extent and projection like 'im_ref'. The mask is used to check if the chosen matching window 

257 position is valid in the sense of useful data. Otherwise, this window position is rejected. 

258 

259 :param CPUs: 

260 number of CPUs to use during calculation of tie point grid (default: None, which means 'all CPUs available') 

261 

262 :param progress: 

263 show progress bars (default: True) 

264 

265 :param v: 

266 verbose mode (default: False) 

267 

268 :param q: 

269 quiet mode (default: False) 

270 

271 :param ignore_errors: 

272 Useful for batch processing. (default: False) 

273 """ 

274 # assertions / input validation 

275 assert gdal.GetDriverByName(fmt_out), "'%s' is not a supported GDAL driver." % fmt_out 

276 if match_gsd and out_gsd: 

277 warnings.warn("'-out_gsd' is ignored because '-match_gsd' is set.\n") 

278 if out_gsd: 

279 assert isinstance(out_gsd, list) and len(out_gsd) == 2, 'out_gsd must be a list with two values.' 

280 

281 self.params = dict([x for x in locals().items() if x[0] != "self" and not x[0].startswith('__')]) 

282 

283 # NOTE: self.imref and self.im2shift are handled completely independent from self.COREG_obj.ref and 

284 # self.COREG_obj.shift. self.COREG_obj.ref and self.COREG_obj.shift are used for shift calculation and 

285 # correction is applied to self.im2shift. 

286 self.imref = GeoArray(im_ref, nodata=nodata[0], progress=progress, q=q) 

287 self.im2shift = GeoArray(im_tgt, nodata=nodata[1], progress=progress, q=q) 

288 self.path_out = path_out # updated by self.set_outpathes 

289 self.fmt_out = fmt_out 

290 self.out_creaOpt = out_crea_options 

291 self._projectDir = projectDir 

292 self.grid_res = grid_res 

293 self.max_points = max_points 

294 self.window_size = window_size 

295 self.max_shift = max_shift 

296 self.max_iter = max_iter 

297 self.tieP_filter_level = tieP_filter_level 

298 self.min_reliability = min_reliability 

299 self.rs_max_outlier = rs_max_outlier 

300 self.rs_tolerance = rs_tolerance 

301 self.rs_random_state = rs_random_state 

302 self.align_grids = align_grids 

303 self.match_gsd = match_gsd 

304 self.out_gsd = out_gsd 

305 self.target_xyGrid = target_xyGrid 

306 self.rspAlg_DS = resamp_alg_deshift # TODO convert integers to strings 

307 self.rspAlg_calc = resamp_alg_calc 

308 self.calc_corners = calc_corners 

309 self.nodata = nodata 

310 self.outFillVal = outFillVal 

311 self.bin_ws = binary_ws 

312 self.force_quadratic_win = force_quadratic_win 

313 self.CPUs = CPUs if CPUs and CPUs <= cpu_count() else cpu_count() 

314 self.path_verbose_out = '' # TODO 

315 self.v = v 

316 self.q = q if not v else False # overridden by v 

317 self.progress = progress if not q else False # overridden by v 

318 self.ignErr = ignore_errors # FIXME this is not yet implemented for COREG_LOCAL 

319 

320 assert self.tieP_filter_level in range(4), 'Invalid tie point filter level.' 

321 assert isinstance(self.imref, GeoArray) and isinstance(self.im2shift, GeoArray), \ 

322 'Something went wrong with the creation of GeoArray instances for reference or target image. The created ' \ 

323 'instances do not seem to belong to the GeoArray class. If you are working in Jupyter Notebook, reset ' \ 

324 'the kernel and try again.' 

325 

326 COREG.__dict__['_set_outpathes'](self, self.imref, self.im2shift) 

327 # make sure that the output directory of coregistered image is the project directory if a project directory is 

328 # given 

329 if path_out and projectDir and os.path.basename(self.path_out): 

330 self.path_out = os.path.join(self.projectDir, os.path.basename(self.path_out)) 

331 

332 # resample input data in case there is a metadata rotation (not handled by AROSICS) 

333 self._check_and_handle_metaRotation() 

334 

335 try: 

336 # ignore_errors must be False because in case COREG init fails, coregistration for the whole scene fails 

337 self.COREG_obj = COREG(self.imref, self.im2shift, 

338 ws=window_size, 

339 footprint_poly_ref=footprint_poly_ref, 

340 footprint_poly_tgt=footprint_poly_tgt, 

341 data_corners_ref=data_corners_ref, 

342 data_corners_tgt=data_corners_tgt, 

343 resamp_alg_calc=self.rspAlg_calc, 

344 calc_corners=calc_corners, 

345 r_b4match=r_b4match, 

346 s_b4match=s_b4match, 

347 max_iter=max_iter, 

348 max_shift=max_shift, 

349 nodata=nodata, 

350 mask_baddata_ref=None, # see below 

351 mask_baddata_tgt=None, 

352 CPUs=self.CPUs, 

353 force_quadratic_win=self.force_quadratic_win, 

354 binary_ws=self.bin_ws, 

355 progress=self.progress, 

356 v=v, 

357 q=q, 

358 ignore_errors=False) 

359 except Exception: 

360 warnings.warn('\nFirst attempt to check the functionality of co-registration failed. Check your ' 

361 'input data and parameters. The following error occurred:', stacklevel=3) 

362 raise 

363 

364 if pyfftw: 

365 self.check_if_fftw_works() 

366 

367 # add bad data mask 

368 # (mask is not added during initialization of COREG object in order to avoid bad data area errors there) 

369 if mask_baddata_ref is not None: 

370 self.COREG_obj.ref.mask_baddata = mask_baddata_ref 

371 if mask_baddata_tgt is not None: 

372 self.COREG_obj.shift.mask_baddata = mask_baddata_tgt 

373 

374 self._tiepoint_grid = None # set by self.tiepoint_grid 

375 self._CoRegPoints_table = None # set by self.CoRegPoints_table 

376 self._coreg_info = None # set by self.coreg_info 

377 self.deshift_results = None # set by self.correct_shifts() 

378 self._success = None # set by self.success property 

379 

380 def _check_and_handle_metaRotation(self): 

381 """Check if the provided input data have a metadata rotation and if yes, correct it AND equalize grids. 

382 

383 In case there is a rotation, the GDAL GeoTransform is not 0 at positions 2 or 4. So far, AROSICS does not 

384 handle such rotations, so the resampling is needed to make things work. The pixel grid equalization is also 

385 done here to avoid a double-resampling (grid would be equalized by COREG.equalize_pixGrids() otherwise). 

386 """ 

387 grid2use = 'ref' if self.im2shift.xgsd <= self.imref.xgsd else 'shift' 

388 

389 if has_metaRotation(self.imref) or has_metaRotation(self.im2shift): 

390 

391 msg = 'The %s image needs to be resampled because it has a row/column rotation in '\ 

392 'its map info which is not handled by AROSICS.' 

393 

394 if grid2use == 'ref': 

395 if has_metaRotation(self.imref): 

396 warnings.warn(msg % 'reference') 

397 self.imref = remove_metaRotation(self.imref) 

398 

399 # resample target to reference image 

400 if not self.q: 

401 print('Adapting the target image pixel grid to the one of the reference image for shift detection.') 

402 self.im2shift.reproject_to_new_grid(prototype=self.imref, CPUs=self.CPUs) 

403 

404 else: 

405 # remove any metadata rotation (a rotation that only exists in the map info) 

406 if has_metaRotation(self.im2shift): 

407 warnings.warn(msg % 'target') 

408 self.im2shift = remove_metaRotation(self.im2shift) 

409 

410 # resample reference to target image 

411 print('Adapting the reference image pixel grid to the one of the target image for shift detection.') 

412 self.imref.reproject_to_new_grid(prototype=self.im2shift, CPUs=self.CPUs) 

413 

414 def check_if_fftw_works(self) -> None: 

415 """Assign the attribute 'fftw_works' to self.COREG_obj by executing shift calculation once with muted output.""" 

416 # calculate global shift once in order to check is fftw works 

417 try: 

418 self.COREG_obj.q = True 

419 self.COREG_obj.v = False 

420 self.COREG_obj.calculate_spatial_shifts() 

421 except RuntimeError: 

422 if self.COREG_obj.fftw_works is not None: 

423 pass 

424 else: 

425 warnings.warn('\nFirst attempt to check if functionality of co-registration failed. Check your ' 

426 'input data and parameters. The following error occurred:', stacklevel=3) 

427 raise 

428 

429 self.COREG_obj.q = self.q 

430 self.COREG_obj.v = self.v 

431 

432 @property 

433 def projectDir(self) -> str: 

434 if self._projectDir: 

435 if len(os.path.split(self._projectDir)) == 1: 

436 return os.path.abspath(os.path.join(os.path.curdir, self._projectDir)) 

437 else: 

438 return os.path.abspath(self._projectDir) 

439 else: 

440 # return a project name that not already has a corresponding folder on disk 

441 root_dir = os.path.dirname(self.im2shift.filePath) if self.im2shift.filePath else os.path.curdir 

442 fold_name = 'UntitledProject_1' 

443 

444 while os.path.isdir(os.path.join(root_dir, fold_name)): 

445 fold_name = '%s_%s' % (fold_name.split('_')[0], int(fold_name.split('_')[-1]) + 1) 

446 

447 self._projectDir = os.path.join(root_dir, fold_name) 

448 return self._projectDir 

449 

450 @property 

451 def tiepoint_grid(self) -> Tie_Point_Grid: 

452 if self._tiepoint_grid: 

453 return self._tiepoint_grid 

454 else: 

455 self.calculate_spatial_shifts() 

456 return self._tiepoint_grid 

457 

458 @property 

459 def CoRegPoints_table(self) -> GeoDataFrame: 

460 """Return a GeoDataFrame containing all the results from coregistration for all points in the tie point grid. 

461 

462 Columns of the GeoDataFrame: 'geometry','POINT_ID','X_IM','Y_IM','X_MAP','Y_MAP','X_WIN_SIZE', 'Y_WIN_SIZE', 

463 'X_SHIFT_PX','Y_SHIFT_PX', 'X_SHIFT_M', 'Y_SHIFT_M', 'ABS_SHIFT' and 'ANGLE' 

464 """ 

465 return self.tiepoint_grid.CoRegPoints_table 

466 

467 @property 

468 def success(self) -> bool: 

469 self._success = self.tiepoint_grid.GCPList != [] 

470 return self._success 

471 

472 def calculate_spatial_shifts(self) -> None: 

473 self._tiepoint_grid = \ 

474 Tie_Point_Grid(self.COREG_obj, self.grid_res, 

475 max_points=self.max_points, 

476 outFillVal=self.outFillVal, 

477 resamp_alg_calc=self.rspAlg_calc, 

478 tieP_filter_level=self.tieP_filter_level, 

479 outlDetect_settings=dict( 

480 min_reliability=self.min_reliability, 

481 rs_max_outlier=self.rs_max_outlier, 

482 rs_tolerance=self.rs_tolerance, 

483 rs_random_state=self.rs_random_state), 

484 dir_out=self.projectDir, 

485 CPUs=self.CPUs, 

486 progress=self.progress, 

487 v=self.v, 

488 q=self.q) 

489 self._tiepoint_grid.get_CoRegPoints_table() 

490 

491 if self.v: 

492 print('Visualizing CoReg points grid...') 

493 self.view_CoRegPoints(figsize=(10, 10)) 

494 

495 def show_image_footprints(self): 

496 """Show a web map containing the calculated footprints and overlap area of the input images. 

497 

498 NOTE: This method is intended to be called from Jupyter Notebook. 

499 """ 

500 return self.COREG_obj.show_image_footprints() 

501 

502 def view_CoRegPoints(self, 

503 shapes2plot: str = 'points', 

504 attribute2plot: str = 'ABS_SHIFT', 

505 cmap: plt.cm = None, 

506 exclude_fillVals: bool = True, 

507 backgroundIm: str = 'tgt', 

508 hide_filtered: bool = True, 

509 figsize: tuple = None, 

510 figsize_multiplier: float = 1, 

511 title: str = '', 

512 vector_scale: float = 1., 

513 savefigPath: str = '', 

514 savefigDPI: int = 96, 

515 showFig: bool = True, 

516 vmin: float = None, 

517 vmax: float = None, 

518 return_map: bool = False 

519 ) -> Optional[Tuple]: 

520 """ 

521 Show a map of the calculated tie point grid with the target image as background. 

522 

523 :param shapes2plot: 'points': plot points representing values of 'attribute2plot' onto the map 

524 'vectors': plot shift vectors onto the map 

525 :param attribute2plot: the attribute of the tie point grid to be shown (default: 'ABS_SHIFT') 

526 :param cmap: a custom color map to be applied to the plotted grid points (default: 'RdYlGn_r') 

527 :param exclude_fillVals: whether to exclude those points of the grid where spatial shift detection failed 

528 :param backgroundIm: whether to use the target or the reference image as map background. Possible 

529 options are 'ref' and 'tgt' (default: 'tgt') 

530 :param hide_filtered: hide all points that have been filtered out according to tie point filter level 

531 :param figsize: size of the figure to be viewed, e.g. (10, 10); automatically estimated if not given 

532 :param figsize_multiplier: if given, the figure size is multiplied with this value 

533 :param title: plot title 

534 :param vector_scale: scale factor for shift vector length (default: 1 -> no scaling) 

535 :param savefigPath: path where to save the figure 

536 :param savefigDPI: DPI resolution of the output figure when saved to disk (default: 96) 

537 :param showFig: whether to show or to hide the figure (default: True) 

538 :param vmin: minimum value of 'attribute2plot' to be included in the figure 

539 :param vmax: maximum value of 'attribute2plot' to be included in the figure 

540 :param return_map: whether to return the figure and axis objects (default: False) 

541 :return: tuple of figure and axis objects or None in case return_map is set to False 

542 """ 

543 from matplotlib import pyplot as plt # noqa 

544 from matplotlib.offsetbox import AnchoredText 

545 from cartopy.crs import PlateCarree 

546 from mpl_toolkits.axes_grid1 import make_axes_locatable 

547 

548 # get background image (reference or target image) 

549 if backgroundIm not in ['tgt', 'ref']: 

550 raise ValueError('backgroundIm') 

551 backgroundIm = self.im2shift if backgroundIm == 'tgt' else self.imref 

552 

553 # make sure the output figure has a reasonable size, also if figsize is not given 

554 if not figsize: 

555 r, c = backgroundIm.shape[:2] 

556 figsize = (8, r / c * 8) if r > c else (c / r * 8, 8) 

557 

558 # apply figsize multiplier 

559 if figsize_multiplier: 

560 if figsize_multiplier < 0: 

561 raise ValueError(figsize_multiplier, 'The figure size multiplier must be a positive finite number.') 

562 figsize = (figsize[0] * figsize_multiplier, figsize[1] * figsize_multiplier) 

563 

564 # get a map showing the background image 

565 fig, ax = backgroundIm.show_map(figsize=figsize, 

566 nodataVal=self.nodata[1], 

567 return_map=True, 

568 band=self.COREG_obj.shift.band4match) 

569 

570 # set figure title 

571 dict_attr_title = dict( 

572 X_WIN_SIZE='size of the matching window in x-direction [pixels]', 

573 Y_WIN_SIZE='size of the matching window in y-direction [pixels]', 

574 X_SHIFT_PX='absolute shifts in x-direction [pixels]', 

575 Y_SHIFT_PX='absolute shifts in y-direction [pixels]', 

576 X_SHIFT_M='absolute shifts in x-direction [map units]', 

577 Y_SHIFT_M='absolute shifts in y-direction [map units]', 

578 ABS_SHIFT='absolute shift vector length [map units]', 

579 ANGLE='shift vector direction [angle in degrees]', 

580 SSIM_BEFORE='structural similarity index before co-registration', 

581 SSIM_AFTER='structural similarity index after co-registration', 

582 SSIM_IMPROVED='structural similarity index improvement through co-registration [yes/no]', 

583 RELIABILITY='reliability of the computed shift vector' 

584 ) 

585 if title: 

586 ax.set_title(title) 

587 elif attribute2plot in dict_attr_title: 

588 ax.set_title(dict_attr_title[attribute2plot], pad=20) 

589 elif attribute2plot in self.CoRegPoints_table.columns: 

590 ax.set_title(attribute2plot) 

591 else: 

592 raise ValueError(attribute2plot, "Invalid value for 'attribute2plot'. Valid values are: %s." 

593 % ", ".join(self.CoRegPoints_table.columns)) 

594 

595 if not self.CoRegPoints_table.empty: 

596 # get GeoDataFrame containing everything needed for plotting 

597 outlierCols = [c for c in self.CoRegPoints_table.columns if 'OUTLIER' in c] 

598 attr2include = ['geometry', attribute2plot] + outlierCols + ['X_SHIFT_M', 'Y_SHIFT_M'] 

599 GDF = self.CoRegPoints_table.loc[self.CoRegPoints_table.X_SHIFT_M != self.outFillVal, attr2include].copy()\ 

600 if exclude_fillVals else self.CoRegPoints_table.loc[:, attr2include] 

601 

602 # get LonLat coordinates for all points 

603 XY = np.array([geom.coords.xy for geom in GDF.geometry]).reshape(-1, 2) 

604 lon, lat = transform_any_prj(self.im2shift.projection, 4326, XY[:, 0], XY[:, 1]) 

605 GDF['Lon'], GDF['Lat'] = lon, lat 

606 

607 # get colors for all points 

608 palette = cmap if cmap is not None else plt.colormaps.get_cmap('RdYlGn_r') 

609 if cmap is None and attribute2plot == 'ANGLE': 

610 import cmocean 

611 palette = getattr(cmocean.cm, 'delta') 

612 

613 if hide_filtered: 

614 if self.tieP_filter_level > 0: 

615 GDF = GDF[GDF.L1_OUTLIER.__eq__(False)].copy() 

616 if self.tieP_filter_level > 1: 

617 GDF = GDF[GDF.L2_OUTLIER.__eq__(False)].copy() 

618 if self.tieP_filter_level > 2: 

619 GDF = GDF[GDF.L3_OUTLIER.__eq__(False)].copy() 

620 else: 

621 marker = 'o' if len(GDF) < 10000 else '.' 

622 common_kw = dict(marker=marker, alpha=1.0, transform=PlateCarree()) 

623 if self.tieP_filter_level > 0: 

624 # flag level 1 outliers 

625 GDF_filt = GDF[GDF.L1_OUTLIER.__eq__(True)].copy() 

626 ax.scatter(GDF_filt['Lon'], GDF_filt['Lat'], c='b', s=250, label='reliability', 

627 **common_kw) 

628 if self.tieP_filter_level > 1: 

629 # flag level 2 outliers 

630 GDF_filt = GDF[GDF.L2_OUTLIER.__eq__(True)].copy() 

631 ax.scatter(GDF_filt['Lon'], GDF_filt['Lat'], c='r', s=150, label='SSIM', 

632 **common_kw) 

633 if self.tieP_filter_level > 2: 

634 # flag level 3 outliers 

635 GDF_filt = GDF[GDF.L3_OUTLIER.__eq__(True)].copy() 

636 ax.scatter(GDF_filt['Lon'], GDF_filt['Lat'], c='y', s=250, label='RANSAC', 

637 **common_kw) 

638 

639 if self.tieP_filter_level > 0: 

640 ax.legend(loc=0, scatterpoints=1) 

641 

642 # plot all points or vectors on top 

643 if not GDF.empty: 

644 vmin_auto, vmax_auto = \ 

645 (np.percentile(GDF[attribute2plot], 0), 

646 np.percentile(GDF[attribute2plot], 98)) \ 

647 if attribute2plot != 'ANGLE' else (0, 360) 

648 vmin = vmin if vmin is not None else vmin_auto 

649 vmax = vmax if vmax is not None else vmax_auto 

650 

651 if shapes2plot == 'vectors': 

652 # plot shift vectors 

653 # doc: https://matplotlib.org/devdocs/api/_as_gen/matplotlib.axes.Axes.quiver.html 

654 mappable = ax.quiver( 

655 GDF['Lon'].values, GDF['Lat'].values, 

656 -GDF['X_SHIFT_M'].values, 

657 -GDF['Y_SHIFT_M'].values, # invert absolute shifts to make arrows point to tgt 

658 GDF[attribute2plot].clip(vmin, vmax), # sets the colors 

659 scale=1200 / vector_scale, # larger values decrease the arrow length 

660 width=.0015, # arrow width (in relation to plot width) 

661 # linewidth=1, # maybe use this to mark outliers instead of scatter points 

662 cmap=palette, 

663 pivot='middle', # position the middle point of the arrows onto the tie point location 

664 transform=PlateCarree() 

665 ) 

666 

667 elif shapes2plot == 'points': 

668 # plot tie points 

669 mappable = ax.scatter( 

670 GDF['Lon'], GDF['Lat'], 

671 c=GDF[attribute2plot], 

672 lw=0, 

673 cmap=palette, 

674 marker='o' if len(GDF) < 10000 else '.', 

675 s=50, 

676 alpha=1.0, 

677 vmin=vmin, 

678 vmax=vmax, 

679 transform=PlateCarree()) 

680 pass 

681 else: 

682 raise ValueError("The parameter 'shapes2plot' must be set to 'vectors' or 'points'. " 

683 "Received %s." % shapes2plot) 

684 

685 # add colorbar 

686 divider = make_axes_locatable(ax) 

687 cax = divider.new_vertical(size="2%", pad=0.4, pack_start=True, 

688 axes_class=plt.Axes # needed because ax is a GeoAxis instance 

689 ) 

690 fig.add_axes(cax) 

691 fig.colorbar(mappable, cax=cax, orientation="horizontal") 

692 

693 # hack to enlarge the figure on the top to avoid cutting off the title (everthing else has no effect) 

694 divider.new_vertical(size="2%", pad=0.4, pack_start=False, axes_class=plt.Axes) 

695 

696 else: 

697 msg = "The map does not contain any tie points \n" \ 

698 "because all the found tie points were flagged as false-positives." 

699 ax.add_artist(AnchoredText(msg, loc='lower center', prop=dict(c='r'))) 

700 

701 if not self.q: 

702 warnings.warn(msg) 

703 

704 else: 

705 msg = "The map does not contain any tie points because no tie points were found at all." 

706 ax.add_artist(AnchoredText(msg, loc='lower center', prop=dict(c='r'))) 

707 

708 if not self.q: 

709 warnings.warn(msg) 

710 

711 # remove white space around the figure 

712 plt.subplots_adjust(top=.95, bottom=.05, right=.95, left=.05) 

713 

714 if savefigPath: 

715 fig.savefig(savefigPath, dpi=savefigDPI, pad_inches=0.1, bbox_inches='tight') 

716 

717 if return_map: 

718 return fig, ax 

719 

720 if showFig and not self.q: 

721 plt.show(block=True) 

722 else: 

723 plt.close(fig) 

724 

725 def view_CoRegPoints_folium(self, attribute2plot: str = 'ABS_SHIFT'): 

726 warnings.warn(UserWarning('This function is still under construction and may not work as expected!')) 

727 assert self.CoRegPoints_table is not None, 'Calculate tie point grid first!' 

728 

729 import folium 

730 import geojson 

731 from folium.raster_layers import ImageOverlay 

732 

733 lon_min, lat_min, lon_max, lat_max = \ 

734 reproject_shapelyGeometry(self.im2shift.box.mapPoly, self.im2shift.projection, 4326).bounds 

735 center_lon, center_lat = (lon_min + lon_max) / 2, (lat_min + lat_max) / 2 

736 

737 # get image to plot 

738 image2plot = self.im2shift[:, :, 0] # FIXME hardcoded band 

739 

740 from py_tools_ds.geo.raster.reproject import warp_ndarray 

741 image2plot, gt, prj = \ 

742 warp_ndarray(image2plot, self.im2shift.geotransform, self.im2shift.projection, 

743 in_nodata=self.nodata[1], out_nodata=self.nodata[1], out_XYdims=(1000, 1000), q=True, 

744 out_prj='epsg:3857') # image must be transformed into web mercator projection 

745 

746 # create map 

747 map_osm = folium.Map(location=[center_lat, center_lon]) # ,zoom_start=3) 

748 # import matplotlib 

749 ImageOverlay( 

750 colormap=lambda x: (1, 0, 0, x), # TODO a colormap must be given 

751 # colormap=matplotlib.cm.gray, # does not work 

752 image=image2plot, bounds=[[lat_min, lon_min], [lat_max, lon_max]], 

753 ).add_to(map_osm) 

754 

755 points_values = self.CoRegPoints_table[['geometry', attribute2plot]] 

756 points_values.geometry.crs = points_values.crs 

757 folium.GeoJson(points_values).add_to(map_osm) 

758 

759 # add overlap polygon 

760 overlapPoly = reproject_shapelyGeometry(self.COREG_obj.overlap_poly, self.im2shift.prj, 4326) 

761 gjs = geojson.Feature(geometry=overlapPoly, properties={}) 

762 folium.GeoJson(gjs).add_to(map_osm) 

763 

764 return map_osm 

765 

766 def _get_updated_map_info_meanShifts(self) -> list: 

767 """Return the updated map info of the target image, shifted on the basis of the mean X/Y shifts.""" 

768 original_map_info = geotransform2mapinfo(self.im2shift.gt, self.im2shift.prj) 

769 updated_map_info = copy(original_map_info) 

770 updated_map_info[3] = str(float(original_map_info[3]) + self.tiepoint_grid.mean_x_shift_map) 

771 updated_map_info[4] = str(float(original_map_info[4]) + self.tiepoint_grid.mean_y_shift_map) 

772 return updated_map_info 

773 

774 @property 

775 def coreg_info(self) -> dict: 

776 """Return a dictionary containing everthing to correct the detected local displacements of the target image.""" 

777 if self._coreg_info: 

778 return self._coreg_info 

779 else: 

780 if not self._tiepoint_grid: 

781 self.calculate_spatial_shifts() 

782 

783 TPG = self._tiepoint_grid 

784 

785 self._coreg_info = { 

786 'GCPList': TPG.GCPList, 

787 'mean_shifts_px': {'x': TPG.mean_x_shift_px if TPG.GCPList else None, 

788 'y': TPG.mean_y_shift_px if TPG.GCPList else None}, 

789 'mean_shifts_map': {'x': TPG.mean_x_shift_map if TPG.GCPList else None, 

790 'y': TPG.mean_y_shift_map if TPG.GCPList else None}, 

791 'updated map info means': self._get_updated_map_info_meanShifts() if TPG.GCPList else None, 

792 'original map info': geotransform2mapinfo(self.imref.gt, self.imref.prj), 

793 'reference projection': self.imref.prj, 

794 'reference geotransform': self.imref.gt, 

795 'reference grid': [[self.imref.gt[0], self.imref.gt[0] + self.imref.gt[1]], 

796 [self.imref.gt[3], self.imref.gt[3] + self.imref.gt[5]]], 

797 'reference extent': {'cols': self.imref.xgsd, 'rows': self.imref.ygsd}, # FIXME not needed anymore 

798 'success': self.success 

799 } 

800 return self.coreg_info 

801 

802 def correct_shifts(self, 

803 max_GCP_count: int = None, 

804 cliptoextent: bool = False, 

805 min_points_local_corr: int = 5 

806 ) -> OrderedDict: 

807 """Perform a local shift correction using all points from the previously calculated tie point grid. 

808 

809 NOTE: Only valid matches are used as GCP points. 

810 

811 :param max_GCP_count: maximum number of GCPs to use 

812 :param cliptoextent: whether to clip the output image to its real extent 

813 :param min_points_local_corr: number of valid tie points, below which a global shift correction is 

814 performed instead of a local correction (global X/Y shift is then computed as 

815 the mean shift of the remaining points)(default: 5 tie points) 

816 :return: 

817 """ 

818 if not self._tiepoint_grid: 

819 self.calculate_spatial_shifts() 

820 

821 if self.tiepoint_grid.GCPList: 

822 if max_GCP_count: 

823 self.coreg_info['GCPList'] = self.coreg_info['GCPList'][:max_GCP_count] 

824 

825 # make sure the correction is applied to the original target image 

826 im2shift = GeoArray(self.params['im_tgt'], nodata=self.nodata[1], progress=self.progress, q=self.q)