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

2 

3# gms_preprocessing, spatial and spectral homogenization of satellite remote sensing data 

4# 

5# Copyright (C) 2020 Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de) 

6# 

7# This software was developed within the context of the GeoMultiSens project funded 

8# by the German Federal Ministry of Education and Research 

9# (project grant code: 01 IS 14 010 A-C). 

10# 

11# This program is free software: you can redistribute it and/or modify it under 

12# the terms of the GNU General Public License as published by the Free Software 

13# Foundation, either version 3 of the License, or (at your option) any later version. 

14# Please note the following exception: `gms_preprocessing` depends on tqdm, which 

15# is distributed under the Mozilla Public Licence (MPL) v2.0 except for the files 

16# "tqdm/_tqdm.py", "setup.py", "README.rst", "MANIFEST.in" and ".gitignore". 

17# Details can be found here: https://github.com/tqdm/tqdm/blob/master/LICENCE. 

18# 

19# This program is distributed in the hope that it will be useful, but WITHOUT 

20# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 

21# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more 

22# details. 

23# 

24# You should have received a copy of the GNU Lesser General Public License along 

25# with this program. If not, see <http://www.gnu.org/licenses/>. 

26 

27""" 

28Level 1B Processor: 

29 

30Detection of global/local geometric displacements. 

31""" 

32 

33 

34import collections 

35import os 

36import time 

37import warnings 

38from datetime import datetime, timedelta 

39 

40import numpy as np 

41from geopandas import GeoDataFrame 

42from shapely.geometry import box 

43import pytz 

44import traceback 

45from typing import Union, TYPE_CHECKING # noqa F401 # flake8 issue 

46 

47from arosics import COREG, DESHIFTER 

48from geoarray import GeoArray 

49from py_tools_ds.geo.coord_grid import is_coord_grid_equal 

50from py_tools_ds.geo.coord_calc import corner_coord_to_minmax 

51from py_tools_ds.geo.coord_trafo import reproject_shapelyGeometry, transform_any_prj 

52from py_tools_ds.geo.projection import prj_equal, EPSG2WKT, WKT2EPSG 

53from py_tools_ds.geo.vector.topology import get_overlap_polygon 

54 

55from ..options.config import GMS_config as CFG 

56from ..model.gms_object import GMS_object 

57from .L1A_P import L1A_object 

58from ..misc import database_tools as DB_T 

59from ..misc import helper_functions as HLP_F 

60from ..misc import path_generator as PG 

61from ..misc.logging import GMS_logger, close_logger 

62from ..misc.locks import DatabaseLock 

63from ..misc.spatial_index_mediator import SpatialIndexMediator 

64from ..misc.definition_dicts import get_GMS_sensorcode, get_outFillZeroSaturated 

65 

66if TYPE_CHECKING: 

67 from shapely.geometry import Polygon # noqa F401 # flake8 issue 

68 from logging import Logger # noqa F401 # flake8 issue 

69 

70__author__ = 'Daniel Scheffler' 

71 

72 

73class Scene_finder(object): 

74 """Scene_finder class to query the postgreSQL database to find a suitable reference scene for co-registration.""" 

75 

76 def __init__(self, src_boundsLonLat, src_AcqDate, src_prj, src_footprint_poly, sceneID_excluded=None, 

77 min_overlap=20, min_cloudcov=0, max_cloudcov=20, plusminus_days=30, plusminus_years=10, logger=None): 

78 # type: (list, datetime, str, Polygon, int, int, int, int, int, int, Logger) -> None 

79 """Initialize Scene_finder. 

80 

81 :param src_boundsLonLat: 

82 :param src_AcqDate: 

83 :param src_prj: 

84 :param src_footprint_poly: 

85 :param sceneID_excluded: 

86 :param min_overlap: minimum overlap of reference scene in percent 

87 :param min_cloudcov: minimum cloud cover of reference scene in percent 

88 :param max_cloudcov: maximum cloud cover of reference scene in percent 

89 :param plusminus_days: maximum time interval between target and reference scene in days 

90 :param plusminus_years: maximum time interval between target and reference scene in years 

91 """ 

92 self.boundsLonLat = src_boundsLonLat 

93 self.src_AcqDate = src_AcqDate 

94 self.src_prj = src_prj 

95 self.src_footprint_poly = src_footprint_poly 

96 self.sceneID_excluded = sceneID_excluded 

97 self.min_overlap = min_overlap 

98 self.min_cloudcov = min_cloudcov 

99 self.max_cloudcov = max_cloudcov 

100 self.plusminus_days = plusminus_days 

101 self.plusminus_years = plusminus_years 

102 self.logger = logger or GMS_logger('ReferenceSceneFinder') 

103 

104 # get temporal constraints 

105 def add_years(dt, years): return dt.replace(dt.year + years) \ 

106 if not (dt.month == 2 and dt.day == 29) else dt.replace(dt.year + years, 3, 1) 

107 self.timeStart = add_years(self.src_AcqDate, -plusminus_years) 

108 timeEnd = add_years(self.src_AcqDate, +plusminus_years) 

109 timeNow = datetime.utcnow().replace(tzinfo=pytz.UTC) 

110 self.timeEnd = timeEnd if timeEnd <= timeNow else timeNow 

111 

112 self.possib_ref_scenes = None # set by self.spatial_query() 

113 self.GDF_ref_scenes = GeoDataFrame() # set by self.spatial_query() 

114 self.ref_scene = None 

115 

116 def __getstate__(self): 

117 """Defines how the attributes of Scene_finder instances are pickled.""" 

118 close_logger(self.logger) 

119 self.logger = None 

120 

121 return self.__dict__ 

122 

123 def __del__(self): 

124 close_logger(self.logger) 

125 self.logger = None 

126 

127 def spatial_query(self, timeout=5): 

128 """Query the postgreSQL database to find possible reference scenes matching the specified criteria. 

129 

130 :param timeout: maximum query duration allowed (seconds) 

131 """ 

132 SpIM = SpatialIndexMediator(host=CFG.spatial_index_server_host, port=CFG.spatial_index_server_port, 

133 timeout=timeout, retries=10) 

134 with DatabaseLock(allowed_slots=1, logger=self.logger): 

135 self.possib_ref_scenes = SpIM.getFullSceneDataForDataset(envelope=self.boundsLonLat, 

136 timeStart=self.timeStart, 

137 timeEnd=self.timeEnd, 

138 minCloudCover=self.min_cloudcov, 

139 maxCloudCover=self.max_cloudcov, 

140 datasetid=CFG.datasetid_spatial_ref, 

141 refDate=self.src_AcqDate, 

142 maxDaysDelta=self.plusminus_days) 

143 

144 if self.possib_ref_scenes: 

145 # fill GeoDataFrame with possible ref scene parameters 

146 GDF = GeoDataFrame(self.possib_ref_scenes, columns=['object']) 

147 GDF['sceneid'] = list(GDF['object'].map(lambda scene: scene.sceneid)) 

148 GDF['acquisitiondate'] = list(GDF['object'].map(lambda scene: scene.acquisitiondate)) 

149 GDF['cloudcover'] = list(GDF['object'].map(lambda scene: scene.cloudcover)) 

150 GDF['polyLonLat'] = list(GDF['object'].map(lambda scene: scene.polyLonLat)) 

151 

152 def LonLat2UTM(polyLL): 

153 return reproject_shapelyGeometry(polyLL, 4326, self.src_prj) 

154 

155 GDF['polyUTM'] = list(GDF['polyLonLat'].map(LonLat2UTM)) 

156 self.GDF_ref_scenes = GDF 

157 

158 def _collect_refscene_metadata(self): 

159 """Collect some reference scene metadata needed for later filtering.""" 

160 GDF = self.GDF_ref_scenes 

161 

162 # get overlap parameter 

163 def get_OL_prms(poly): return get_overlap_polygon(poly, self.src_footprint_poly) 

164 

165 GDF['overlapParams'] = list(GDF['polyLonLat'].map(get_OL_prms)) 

166 GDF['overlap area'] = list(GDF['overlapParams'].map(lambda OL_prms: OL_prms['overlap area'])) 

167 GDF['overlap percentage'] = list(GDF['overlapParams'].map(lambda OL_prms: OL_prms['overlap percentage'])) 

168 GDF['overlap poly'] = list(GDF['overlapParams'].map(lambda OL_prms: OL_prms['overlap poly'])) 

169 del GDF['overlapParams'] 

170 

171 # get processing level of reference scenes 

172 procL = GeoDataFrame( 

173 DB_T.get_info_from_postgreSQLdb(CFG.conn_database, 'scenes_proc', ['sceneid', 'proc_level'], 

174 {'sceneid': list(GDF.sceneid)}), columns=['sceneid', 'proc_level']) 

175 GDF = GDF.merge(procL, on='sceneid', how='left') 

176 GDF = GDF.where(GDF.notnull(), None) # replace NaN values with None 

177 

178 # get path of binary file 

179 def get_path_binary(GDF_row): 

180 return PG.path_generator(scene_ID=GDF_row['sceneid'], proc_level=GDF_row['proc_level']) \ 

181 .get_path_imagedata() if GDF_row['proc_level'] else None 

182 GDF['path_ref'] = GDF.apply(lambda GDF_row: get_path_binary(GDF_row), axis=1) 

183 GDF['refDs_exists'] = list(GDF['path_ref'].map(lambda p: os.path.exists(p) if p else False)) 

184 

185 # check if a proper entity ID can be gathered from database 

186 eID = GeoDataFrame(DB_T.get_info_from_postgreSQLdb(CFG.conn_database, 'scenes', ['id', 'entityid'], 

187 {'id': list(GDF.sceneid)}), columns=['sceneid', 'entityid']) 

188 GDF = GDF.merge(eID, on='sceneid', how='left') 

189 self.GDF_ref_scenes = GDF.where(GDF.notnull(), None) 

190 

191 def _filter_excluded_sceneID(self): 

192 """Filter reference scene with the same scene ID like the target scene.""" 

193 GDF = self.GDF_ref_scenes 

194 if not GDF.empty: 

195 self.logger.info('Same ID filter: Excluding scene with the same ID like the target scene.') 

196 self.GDF_ref_scenes = GDF.loc[GDF['sceneid'] != self.sceneID_excluded] 

197 self.logger.info('%s scenes => %s scenes' % (len(GDF), len(self.GDF_ref_scenes))) 

198 

199 def _filter_by_overlap(self): 

200 """Filter all scenes with less spatial overlap than self.min_overlap.""" 

201 GDF = self.GDF_ref_scenes 

202 if not GDF.empty: 

203 self.logger.info('Overlap filter: Excluding all scenes with less than %s percent spatial overlap.' 

204 % self.min_overlap) 

205 self.GDF_ref_scenes = GDF.loc[GDF['overlap percentage'] >= self.min_overlap] 

206 self.logger.info('%s scenes => %s scenes' % (len(GDF), len(self.GDF_ref_scenes))) 

207 

208 def _filter_by_proc_status(self): 

209 """Filter all scenes that have not been processed before according to proc. status (at least L1A is needed).""" 

210 GDF = self.GDF_ref_scenes 

211 if not GDF.empty: 

212 self.logger.info('Processing level filter: Exclude all scenes that have not been processed before ' 

213 'according to processing status (at least L1A is needed).') 

214 self.GDF_ref_scenes = GDF[GDF['proc_level'].notnull()] 

215 self.logger.info('%s scenes => %s scenes' % (len(GDF), len(self.GDF_ref_scenes))) 

216 

217 def _filter_by_dataset_existence(self): 

218 """Filter all scenes where no processed data can be found on fileserver.""" 

219 GDF = self.GDF_ref_scenes 

220 if not GDF.empty: 

221 self.logger.info('Existence filter: Excluding all scenes where no processed data have been found.') 

222 self.GDF_ref_scenes = GDF[GDF['refDs_exists']] 

223 self.logger.info('%s scenes => %s scenes' % (len(GDF), len(self.GDF_ref_scenes))) 

224 

225 def _filter_by_entity_ID_availability(self): 

226 """Filter all scenes where no proper entity ID can be found in the database (database errors).""" 

227 GDF = self.GDF_ref_scenes 

228 if not GDF.empty: 

229 self.logger.info('DB validity filter: Exclude all scenes where no proper entity ID can be found in the ' 

230 'database (database errors).') 

231 self.GDF_ref_scenes = GDF[GDF['entityid'].notnull()] 

232 self.logger.info('%s scenes => %s scenes' % (len(GDF), len(self.GDF_ref_scenes))) 

233 

234 def _filter_by_projection(self): 

235 """Filter all scenes that have a different projection than the target image.""" 

236 GDF = self.GDF_ref_scenes[self.GDF_ref_scenes.refDs_exists] 

237 if not GDF.empty: 

238 # compare projections of target and reference image 

239 GDF['prj_equal'] = \ 

240 list(GDF['path_ref'].map(lambda path_ref: prj_equal(self.src_prj, GeoArray(path_ref).prj))) 

241 

242 self.logger.info('Projection filter: Exclude all scenes that have a different projection than the target ' 

243 'image.') 

244 self.GDF_ref_scenes = GDF[GDF['prj_equal']] 

245 self.logger.info('%s scenes => %s scenes' % (len(GDF), len(self.GDF_ref_scenes))) 

246 

247 def choose_ref_scene(self): 

248 """Choose reference scene with minimum cloud cover and maximum overlap.""" 

249 if self.possib_ref_scenes: 

250 # First, collect some relavant reference scene metadata 

251 self._collect_refscene_metadata() 

252 

253 # Filter possible scenes by running all filter functions 

254 self._filter_excluded_sceneID() 

255 self._filter_by_overlap() 

256 self._filter_by_proc_status() 

257 self._filter_by_dataset_existence() 

258 self._filter_by_entity_ID_availability() 

259 self._filter_by_projection() 

260 

261 # Choose the reference scene out of the filtered DataFrame 

262 if not self.GDF_ref_scenes.empty: 

263 GDF = self.GDF_ref_scenes 

264 GDF = GDF[GDF['cloudcover'] == GDF['cloudcover'].min()] 

265 GDF = GDF[GDF['overlap percentage'] == GDF['overlap percentage'].max()] 

266 

267 if not GDF.empty: 

268 GDF_res = GDF.iloc[0] 

269 return ref_Scene(GDF_res) 

270 else: 

271 return None 

272 

273 

274class ref_Scene: 

275 def __init__(self, GDF_record): 

276 self.scene_ID = int(GDF_record['sceneid']) 

277 self.entity_ID = GDF_record['entityid'] 

278 self.AcqDate = GDF_record['acquisitiondate'] 

279 self.cloudcover = GDF_record['cloudcover'] 

280 self.polyLonLat = GDF_record['polyLonLat'] 

281 self.polyUTM = GDF_record['polyUTM'] 

282 self.proc_level = GDF_record['proc_level'] 

283 self.filePath = GDF_record['path_ref'] 

284 

285 

286class L1B_object(L1A_object): 

287 def __init__(self, L1A_obj=None): 

288 

289 super(L1B_object, self).__init__() 

290 

291 # set defaults 

292 self._spatRef_available = None 

293 self.spatRef_scene = None # set by self.get_spatial_reference_scene() 

294 self.deshift_results = collections.OrderedDict() 

295 

296 if L1A_obj: 

297 # populate attributes 

298 [setattr(self, key, value) for key, value in L1A_obj.__dict__.items()] 

299 

300 self.proc_level = 'L1B' 

301 self.proc_status = 'initialized' 

302 

303 @property 

304 def spatRef_available(self): 

305 if self._spatRef_available is not None: 

306 return self._spatRef_available 

307 else: 

308 self.get_spatial_reference_scene() 

309 return self._spatRef_available 

310 

311 @spatRef_available.setter 

312 def spatRef_available(self, spatRef_available): 

313 self._spatRef_available = spatRef_available 

314 

315 def get_spatial_reference_scene(self): 

316 boundsLonLat = corner_coord_to_minmax(self.trueDataCornerLonLat) 

317 footprint_poly = HLP_F.CornerLonLat_to_shapelyPoly(self.trueDataCornerLonLat) 

318 RSF = Scene_finder(boundsLonLat, self.acq_datetime, self.MetaObj.projection, 

319 footprint_poly, self.scene_ID, 

320 min_overlap=CFG.spatial_ref_min_overlap, 

321 min_cloudcov=CFG.spatial_ref_min_cloudcov, 

322 max_cloudcov=CFG.spatial_ref_max_cloudcov, 

323 plusminus_days=CFG.spatial_ref_plusminus_days, 

324 plusminus_years=CFG.spatial_ref_plusminus_years, 

325 logger=self.logger) 

326 

327 # run spatial query 

328 self.logger.info('Querying database in order to find a suitable reference scene for co-registration.') 

329 RSF.spatial_query(timeout=5) 

330 if RSF.possib_ref_scenes: 

331 self.logger.info('Query result: %s reference scenes with matching metadata.' % len(RSF.possib_ref_scenes)) 

332 

333 # try to get a spatial reference scene by applying some filter criteria 

334 self.spatRef_scene = RSF.choose_ref_scene() # type: Union[ref_Scene, None] 

335 if self.spatRef_scene: 

336 self.spatRef_available = True 

337 self.logger.info('Found a suitable reference image for coregistration: scene ID %s (entity ID %s).' 

338 % (self.spatRef_scene.scene_ID, self.spatRef_scene.entity_ID)) 

339 else: 

340 self.spatRef_available = False 

341 self.logger.warning('No scene fulfills all requirements to serve as spatial reference for scene %s ' 

342 '(entity ID %s). Coregistration impossible.' % (self.scene_ID, self.entity_ID)) 

343 

344 else: 

345 self.logger.warning('Spatial query returned no matches. Coregistration impossible.') 

346 self.spatRef_available = False 

347 

348 def _get_reference_image_params_pgSQL(self): 

349 # TODO implement earlier version of this function as a backup for SpatialIndexMediator 

350 """postgreSQL query: get IDs of overlapping scenes and select most suitable scene_ID 

351 (with respect to DGM, cloud cover""" 

352 warnings.warn('_get_reference_image_params_pgSQL is deprecated an will not work anymore.', DeprecationWarning) 

353 

354 # vorab-check anhand wolkenmaske, welche region von im2shift überhaupt für shift-corr tauglich ist 

355 # -> diese region als argument in postgresql abfrage 

356 # scene_ID = 14536400 # LE71510322000093SGS00 im2shift 

357 

358 # set query conditions 

359 min_overlap = 20 # % 

360 max_cloudcov = 20 # % 

361 plusminus_days = 30 

362 AcqDate = self.im2shift_objDict['acquisition_date'] 

363 date_minmax = [AcqDate - timedelta(days=plusminus_days), AcqDate + timedelta(days=plusminus_days)] 

364 dataset_cond = 'datasetid=%s' % CFG.datasetid_spatial_ref 

365 cloudcov_cond = 'cloudcover < %s' % max_cloudcov 

366 # FIXME cloudcover noch nicht für alle scenes im proc_level METADATA verfügbar 

367 dayrange_cond = "(EXTRACT(MONTH FROM scenes.acquisitiondate), EXTRACT(DAY FROM scenes.acquisitiondate)) " \ 

368 "BETWEEN (%s, %s) AND (%s, %s)" \ 

369 % (date_minmax[0].month, date_minmax[0].day, date_minmax[1].month, date_minmax[1].day) 

370 # TODO weitere Kriterien einbauen! 

371 

372 def query_scenes(condlist): 

373 return DB_T.get_overlapping_scenes_from_postgreSQLdb( 

374 CFG.conn_database, 

375 table='scenes', 

376 tgt_corners_lonlat=self.trueDataCornerLonLat, 

377 conditions=condlist, 

378 add_cmds='ORDER BY scenes.cloudcover ASC', 

379 timeout=30000) 

380 conds_descImportance = [dataset_cond, cloudcov_cond, dayrange_cond] 

381 

382 self.logger.info('Querying database in order to find a suitable reference scene for co-registration.') 

383 

384 count, filt_overlap_scenes = 0, [] 

385 while not filt_overlap_scenes: 

386 if count == 0: 

387 # search within already processed scenes 

388 # das ist nur Ergebnis aus scenes_proc 

389 # -> dort liegt nur eine referenz, wenn die szene schon bei CFG.job-Beginn in Datensatzliste drin war 

390 res = DB_T.get_overlapping_scenes_from_postgreSQLdb( 

391 CFG.conn_database, 

392 tgt_corners_lonlat=self.trueDataCornerLonLat, 

393 conditions=['datasetid=%s' % CFG.datasetid_spatial_ref], 

394 add_cmds='ORDER BY scenes.cloudcover ASC', 

395 timeout=25000) 

396 filt_overlap_scenes = self._sceneIDList_to_filt_overlap_scenes([i[0] for i in res[:50]], 20.) 

397 

398 else: 

399 # search within complete scenes table using less filter criteria with each run 

400 # TODO: Daniels Index einbauen, sonst bei wachsender Tabellengröße evtl. irgendwann zu langsam 

401 res = query_scenes(conds_descImportance) 

402 filt_overlap_scenes = self._sceneIDList_to_filt_overlap_scenes([i[0] for i in res[:50]], min_overlap) 

403 

404 if len(conds_descImportance) > 1: # FIXME anderer Referenzsensor? 

405 del conds_descImportance[-1] 

406 else: # reduce min_overlap to 10 percent if there are overlapping scenes 

407 if res: 

408 min_overlap = 10 

409 filt_overlap_scenes = \ 

410 self._sceneIDList_to_filt_overlap_scenes([i[0] for i in res[:50]], min_overlap) 

411 

412 # raise warnings if no match found 

413 if not filt_overlap_scenes: 

414 if not res: 

415 warnings.warn('No reference scene found for %s (scene ID %s). Coregistration of this scene ' 

416 'failed.' % (self.baseN, self.scene_ID)) 

417 else: 

418 warnings.warn('Reference scenes for %s (scene ID %s) have been found but none has more ' 

419 'than %s percent overlap. Coregistration of this scene failed.' 

420 % (self.baseN, self.scene_ID, min_overlap)) 

421 break 

422 count += 1 

423 

424 if filt_overlap_scenes: 

425 ref_available = False 

426 for count, sc in enumerate(filt_overlap_scenes): 

427 if count == 2: # FIXME Abbuch schon bei 3. Szene? 

428 warnings.warn('No reference scene for %s (scene ID %s) available. ' 

429 'Coregistration of this scene failed.' % (self.baseN, self.scene_ID)) 

430 break 

431 

432 # start download of scene data not available and start L1A processing 

433 def dl_cmd(scene_ID): print('%s %s %s' % ( 

434 CFG.java_commands['keyword'].strip(), # FIXME CFG.java_commands is deprecated 

435 CFG.java_commands["value_download"].strip(), scene_ID)) 

436 

437 path = PG.path_generator(scene_ID=sc['scene_ID']).get_path_imagedata() 

438 

439 if not os.path.exists(path): 

440 # add scene 2 download to scenes_jobs.missing_scenes 

441 

442 # print JAVA download command 

443 t_dl_start = time.time() 

444 dl_cmd(sc['scene_ID']) 

445 

446 # check if scene is downloading 

447 download_start_timeout = 5 # seconds 

448 # set timout for external processing 

449 # -> DEPRECATED BECAUSE CREATION OF EXTERNAL CFG WITHIN CFG IS NOT ALLOWED 

450 processing_timeout = 5 # seconds # FIXME increase timeout if processing is really started 

451 proc_level = None 

452 while True: 

453 proc_level_chk = DB_T.get_info_from_postgreSQLdb( 

454 CFG.conn_database, 'scenes', ['proc_level'], {'id': sc['scene_ID']})[0][0] 

455 if proc_level != proc_level_chk: 

456 print('Reference scene %s, current processing level: %s' % (sc['scene_ID'], proc_level_chk)) 

457 proc_level = proc_level_chk 

458 if proc_level_chk in ['SCHEDULED', 'METADATA'] and \ 

459 time.time() - t_dl_start >= download_start_timeout: 

460 warnings.warn('Download of reference scene %s (entity ID %s) timed out. ' 

461 'Coregistration of this scene failed.' % (self.baseN, self.scene_ID)) 

462 break 

463 if proc_level_chk == 'L1A': 

464 ref_available = True 

465 break 

466 elif proc_level_chk == 'DOWNLOADED' and time.time() - t_dl_start >= processing_timeout: 

467 # proc_level='DOWNLOADED' but scene has not been processed 

468 warnings.warn('L1A processing of reference scene %s (entity ID %s) timed out. ' 

469 'Coregistration of this scene failed.' % (self.baseN, self.scene_ID)) 

470 break 

471 # DB_T.set_info_in_postgreSQLdb(CFG.conn_database,'scenes', 

472 # {'proc_level':'METADATA'},{'id':sc['scene_ID']}) 

473 

474 time.sleep(5) 

475 else: 

476 ref_available = True 

477 

478 if not ref_available: 

479 continue 

480 else: 

481 self.path_imref = path 

482 self.imref_scene_ID = sc['scene_ID'] 

483 self.imref_footprint_poly = sc['scene poly'] 

484 self.overlap_poly = sc['overlap poly'] 

485 self.overlap_percentage = sc['overlap percentage'] 

486 self.overlap_area = sc['overlap area'] 

487 

488 query_res = DB_T.get_info_from_postgreSQLdb(CFG.conn_database, 'scenes', ['entityid'], 

489 {'id': self.imref_scene_ID}, records2fetch=1) 

490 assert query_res != [], 'No entity-ID found for scene number %s' % self.imref_scene_ID 

491 self.imref_entity_ID = query_res[0][0] # [('LC81510322013152LGN00',)] 

492 break 

493 self.logger.close() 

494 

495 def _sceneIDList_to_filt_overlap_scenes(self, sceneIDList, min_overlap): 

496 """find reference scenes that cover at least 20% of the scene with the given ID 

497 ONLY FIRST 50 scenes are considered""" 

498 

499 # t0 = time.time() 

500 dict_sceneID_poly = [{'scene_ID': ID, 'scene poly': HLP_F.scene_ID_to_shapelyPolygon(ID)} 

501 for ID in sceneIDList] # always returns LonLot polygons 

502 

503 # get overlap polygons and their parameters. result: [{overlap poly, overlap percentage, overlap area}] 

504 for dic in dict_sceneID_poly: # input: dicts {scene_ID, ref_poly} 

505 dict_overlap_poly_params = get_overlap_polygon(dic['scene poly'], self.arr.footprint_poly) 

506 dic.update(dict_overlap_poly_params) # adds {overlap poly, overlap percentage, overlap area} 

507 # print('polygon creation time', time.time()-t0) 

508 

509 # filter those scene_IDs out where overlap percentage is below 20% 

510 if min_overlap: 

511 filt_overlap_scenes = [scene for scene in dict_sceneID_poly if scene['overlap percentage'] >= min_overlap] 

512 else: 

513 filt_overlap_scenes = dict_sceneID_poly 

514 

515 return filt_overlap_scenes 

516 

517 def get_opt_bands4matching(self, target_cwlPos_nm=550): 

518 """Automatically determines the optimal bands used für fourier shift theorem matching 

519 

520 :param target_cwlPos_nm: the desired wavelength used for matching 

521 """ 

522 # get GMS_object for reference scene 

523 path_gmsFile = PG.path_generator(scene_ID=self.spatRef_scene.scene_ID).get_path_gmsfile() 

524 ref_obj = GMS_object.from_disk((path_gmsFile, ['cube', None])) 

525 

526 # get spectral characteristics 

527 ref_cwl = [float(ref_obj.MetaObj.CWL[bN]) for bN in ref_obj.MetaObj.LayerBandsAssignment] 

528 shift_cwl = [float(self.MetaObj.CWL[bN]) for bN in self.MetaObj.LayerBandsAssignment] 

529 shift_fwhm = [float(self.MetaObj.FWHM[bN]) for bN in self.MetaObj.LayerBandsAssignment] 

530 

531 # exclude cirrus/oxygen band of Landsat-8/Sentinel-2 

532 shift_bbl, ref_bbl = [False] * len(shift_cwl), [False] * len(ref_cwl) # bad band lists 

533 for GMS_obj, s_r, bbl in zip([self, ref_obj], ['shift', 'ref'], [shift_bbl, ref_bbl]): 

534 GMS_obj.GMS_identifier.logger = None # set a dummy value in order to avoid Exception 

535 sensorcode = get_GMS_sensorcode(GMS_obj.GMS_identifier) 

536 if sensorcode in ['LDCM', 'S2A', 'S2B'] and '9' in GMS_obj.LayerBandsAssignment: 

537 bbl[GMS_obj.LayerBandsAssignment.index('9')] = True 

538 if sensorcode in ['S2A', 'S2B'] and '10' in GMS_obj.LayerBandsAssignment: 

539 bbl[GMS_obj.LayerBandsAssignment.index('10')] = True 

540 

541 # cwl_overlap = (max(min(shift_cwl),min(ref_cwl)), min(max(shift_cwl),max(ref_cwl))) # -> (min wvl, max wvl) 

542 # find matching band of reference image for each band of image to be shifted 

543 match_dic = collections.OrderedDict() 

544 for idx, cwl, fwhm in zip(range(len(shift_cwl)), shift_cwl, shift_fwhm): 

545 if shift_bbl[idx]: 

546 continue # skip cwl if it is declared as bad band above 

547 

548 def is_inside(r_cwl, s_cwl, s_fwhm): return s_cwl - s_fwhm / 2 < r_cwl < s_cwl + s_fwhm / 2 

549 

550 matching_r_cwls = [r_cwl for i, r_cwl in enumerate(ref_cwl) if 

551 is_inside(r_cwl, cwl, fwhm) and not ref_bbl[i]] 

552 if matching_r_cwls: 

553 match_dic[cwl] = matching_r_cwls[0] if len(matching_r_cwls) else \ 

554 matching_r_cwls[np.abs(np.array(matching_r_cwls) - cwl).argmin()] 

555 

556 # set bands4 match based on the above results 

557 poss_cwls = [cwl for cwl in shift_cwl if cwl in match_dic] 

558 if poss_cwls: 

559 if not target_cwlPos_nm: 

560 shift_band4match = shift_cwl.index(poss_cwls[0]) + 1 # first possible shift band 

561 ref_band4match = ref_cwl.index(match_dic[poss_cwls[0]]) + 1 # matching reference band 

562 else: # target_cwlPos is given 

563 closestWvl_to_target = poss_cwls[np.abs(np.array(poss_cwls) - target_cwlPos_nm).argmin()] 

564 shift_band4match = shift_cwl.index(closestWvl_to_target) + 1 # the shift band closest to target 

565 ref_band4match = ref_cwl.index(match_dic[closestWvl_to_target]) + 1 # matching ref 

566 else: # all reference bands are outside of shift-cwl +- fwhm/2 

567 self.logger.warning('Optimal bands for matching could not be automatically determined. ' 

568 'Choosing first band of each image.') 

569 shift_band4match = 1 

570 ref_band4match = 1 

571 

572 self.logger.info( 

573 'Target band for matching: %s (%snm)' % (shift_band4match, shift_cwl[shift_band4match - 1])) 

574 self.logger.info( 

575 'Reference band for matching: %s (%snm)' % (ref_band4match, ref_cwl[ref_band4match - 1])) 

576 

577 return ref_band4match, shift_band4match 

578 

579 def compute_global_shifts(self): 

580 spatIdxSrv_status = os.environ['GMS_SPAT_IDX_SRV_STATUS'] if 'GMS_SPAT_IDX_SRV_STATUS' in os.environ else True 

581 

582 if spatIdxSrv_status == 'unavailable': 

583 self.logger.warning('Coregistration skipped due to unavailable Spatial Index Mediator Server!"') 

584 

585 elif CFG.skip_coreg: 

586 self.logger.warning('Coregistration skipped according to user configuration.') 

587 

588 elif self.coreg_needed and self.spatRef_available: 

589 self.coreg_info.update({'reference scene ID': self.spatRef_scene.scene_ID}) 

590 self.coreg_info.update({'reference entity ID': self.spatRef_scene.entity_ID}) 

591 

592 geoArr_ref = GeoArray(self.spatRef_scene.filePath) 

593 geoArr_shift = GeoArray(self.arr) 

594 r_b4match, s_b4match = self.get_opt_bands4matching(target_cwlPos_nm=CFG.coreg_band_wavelength_for_matching) 

595 coreg_kwargs = dict( 

596 r_b4match=r_b4match, 

597 s_b4match=s_b4match, 

598 align_grids=True, # FIXME not needed here 

599 match_gsd=True, # FIXME not needed here 

600 max_shift=CFG.coreg_max_shift_allowed, 

601 ws=CFG.coreg_window_size, 

602 data_corners_ref=[[x, y] for x, y in self.spatRef_scene.polyUTM.convex_hull.exterior.coords], 

603 data_corners_tgt=[transform_any_prj(EPSG2WKT(4326), self.MetaObj.projection, x, y) 

604 for x, y in HLP_F.reorder_CornerLonLat(self.trueDataCornerLonLat)], 

605 nodata=(get_outFillZeroSaturated(geoArr_ref.dtype)[0], 

606 get_outFillZeroSaturated(geoArr_shift.dtype)[0]), 

607 ignore_errors=True, 

608 v=False, 

609 q=True 

610 ) 

611 

612 # initialize COREG object 

613 try: 

614 COREG_obj = COREG(geoArr_ref, geoArr_shift, **coreg_kwargs) 

615 except Exception as e: 

616 COREG_obj = None 

617 self.logger.error('\nAn error occurred during coregistration. BE AWARE THAT THE SCENE %s ' 

618 '(ENTITY ID %s) HAS NOT BEEN COREGISTERED! Error message was: \n%s\n' 

619 % (self.scene_ID, self.entity_ID, repr(e))) 

620 self.logger.error(traceback.format_exc()) 

621 # TODO include that in the job summary 

622 

623 # calculate_spatial_shifts 

624 if COREG_obj: 

625 COREG_obj.calculate_spatial_shifts() 

626 

627 self.coreg_info.update( 

628 COREG_obj.coreg_info) # no clipping to trueCornerLonLat until here -> only shift correction 

629 self.coreg_info.update({'shift_reliability': COREG_obj.shift_reliability}) 

630 

631 if COREG_obj.success: 

632 self.coreg_info['success'] = True 

633 self.logger.info("Calculated map shifts (X,Y): %s / %s" 

634 % (COREG_obj.x_shift_map, 

635 COREG_obj.y_shift_map)) # FIXME direkt in calculate_spatial_shifts loggen 

636 self.logger.info("Reliability of calculated shift: %.1f percent" % COREG_obj.shift_reliability) 

637 

638 else: 

639 # TODO add database entry with error hint 

640 [self.logger.error('ERROR during coregistration of scene %s (entity ID %s):\n%s' 

641 % (self.scene_ID, self.entity_ID, err)) for err in COREG_obj.tracked_errors] 

642 

643 else: 

644 if self.coreg_needed: 

645 self.logger.warning('Coregistration skipped because no suitable reference scene is available or ' 

646 'spatial query failed.') 

647 else: 

648 self.logger.info('Coregistration of scene %s (entity ID %s) skipped because target dataset ID equals ' 

649 'reference dataset ID.' % (self.scene_ID, self.entity_ID)) 

650 

651 def correct_spatial_shifts(self, cliptoextent=True, clipextent=None, clipextent_prj=None, v=False): 

652 # type: (bool, list, any, bool) -> None 

653 """Corrects the spatial shifts calculated by self.compute_global_shifts(). 

654 

655 :param cliptoextent: whether to clip the output to the given extent 

656 :param clipextent: list of XY-coordinate tuples giving the target extent (if not given and cliptoextent is 

657 True, the 'trueDataCornerLonLat' attribute of the GMS object is used 

658 :param clipextent_prj: WKT projection string or EPSG code of the projection for the coordinates in clipextent 

659 :param v: 

660 :return: 

661 """ 

662 

663 # cliptoextent is automatically True if an extent is given 

664 cliptoextent = cliptoextent if not clipextent else True 

665 

666 if cliptoextent or self.resamp_needed or (self.coreg_needed and self.coreg_info['success']): 

667 

668 # get target bounds # TODO implement boxObj call instead here 

669 if not clipextent: 

670 trueDataCornerUTM = [transform_any_prj(EPSG2WKT(4326), self.MetaObj.projection, x, y) 

671 for x, y in self.trueDataCornerLonLat] 

672 xmin, xmax, ymin, ymax = corner_coord_to_minmax(trueDataCornerUTM) 

673 mapBounds = box(xmin, ymin, xmax, ymax).bounds 

674 else: 

675 assert clipextent_prj, \ 

676 "'clipextent_prj' must be given together with 'clipextent'. Received only 'clipextent'." 

677 clipextent_UTM = clipextent if prj_equal(self.MetaObj.projection, clipextent_prj) else \ 

678 [transform_any_prj(clipextent_prj, self.MetaObj.projection, x, y) for x, y in clipextent] 

679 xmin, xmax, ymin, ymax = corner_coord_to_minmax(clipextent_UTM) 

680 mapBounds = box(xmin, ymin, xmax, ymax).bounds 

681 

682 # correct shifts and clip to extent 

683 # ensure self.masks exists (does not exist in case of inmem_serialization mode because 

684 # then self.fill_from_disk() is skipped) 

685 if not hasattr(self, 'masks') or self.masks is None: 

686 self.build_combined_masks_array() # creates self.masks and self.masks_meta 

687 

688 # exclude self.mask_nodata, self.mask_clouds from warping 

689 del self.mask_nodata, self.mask_clouds 

690 

691 attributes2deshift = [attrname for attrname in 

692 ['arr', 'masks', 'dem', 'ac_errors', 'mask_clouds_confidence'] 

693 if getattr(self, '_%s' % attrname) is not None] 

694 for attrname in attributes2deshift: 

695 geoArr = getattr(self, attrname) 

696 

697 # do some logging 

698 if self.coreg_needed: 

699 if self.coreg_info['success']: 

700 self.logger.info("Correcting spatial shifts for attribute '%s'..." % attrname) 

701 elif cliptoextent and is_coord_grid_equal( 

702 geoArr.gt, CFG.spatial_ref_gridx, CFG.spatial_ref_gridy): 

703 self.logger.info("Attribute '%s' has only been clipped to it's extent because no valid " 

704 "shifts have been detected and the pixel grid equals the target grid." 

705 % attrname) 

706 elif self.resamp_needed: 

707 self.logger.info("Resampling attribute '%s' to target grid..." % attrname) 

708 elif self.resamp_needed: 

709 self.logger.info("Resampling attribute '%s' to target grid..." % attrname) 

710 

711 # correct shifts 

712 DS = DESHIFTER(geoArr, self.coreg_info, 

713 target_xyGrid=[CFG.spatial_ref_gridx, CFG.spatial_ref_gridy], 

714 cliptoextent=cliptoextent, 

715 clipextent=mapBounds, 

716 align_grids=True, 

717 resamp_alg='nearest' if attrname == 'masks' else CFG.spatial_resamp_alg, 

718 CPUs=None if CFG.allow_subMultiprocessing else 1, 

719 progress=True if v else False, 

720 q=True, 

721 v=v) 

722 DS.correct_shifts() 

723 

724 # update coreg_info 

725 if attrname == 'arr': 

726 self.coreg_info['is shifted'] = DS.is_shifted 

727 self.coreg_info['is resampled'] = DS.is_resampled 

728 

729 # update geoinformations and array shape related attributes 

730 self.logger.info("Updating geoinformations of '%s' attribute..." % attrname) 

731 if attrname == 'arr': 

732 self.MetaObj.map_info = DS.updated_map_info 

733 self.MetaObj.projection = EPSG2WKT(WKT2EPSG(DS.updated_projection)) 

734 self.shape_fullArr = DS.arr_shifted.shape 

735 self.MetaObj.rows, self.MetaObj.cols = DS.arr_shifted.shape[:2] 

736 else: 

737 self.masks_meta['map info'] = DS.updated_map_info 

738 self.masks_meta['coordinate system string'] = EPSG2WKT(WKT2EPSG(DS.updated_projection)) 

739 self.masks_meta['lines'], self.masks_meta['samples'] = DS.arr_shifted.shape[:2] 

740 

741 # NOTE: mask_nodata and mask_clouds are updated later by L2A_map mapper function (module pipeline) 

742 

743 # update the GeoArray instance without loosing its inherent metadata (nodata, ...) 

744 geoArr.arr, geoArr.gt = DS.GeoArray_shifted.arr, DS.GeoArray_shifted.gt 

745 if not geoArr.prj: 

746 geoArr.prj = DS.GeoArray_shifted.prj 

747 # setattr(self,attrname, DS.GeoArray_shifted) # NOTE: don't set array earlier because setter will also 

748 # # update arr.gt/.prj/.nodata from MetaObj 

749 

750 self.resamp_needed = False 

751 self.coreg_needed = False 

752 

753 # recreate self.masks_nodata and self.mask_clouds from self.masks 

754 self.mask_nodata = self.mask_nodata 

755 self.mask_clouds = self.mask_clouds 

756 # FIXME move functionality of self.masks only to output writer and remove self.masks completely