Coverage for gms_preprocessing/algorithms/L1C_P.py : 82%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# -*- coding: utf-8 -*-
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/>.
27"""Level 1C Processor: Atmospheric correction of TOA-reflectance data."""
29import warnings
30import re
31import logging
32import dill
33import traceback
34from typing import List # noqa F401 # flake8 issue
35from time import time
36import os
37import timeout_decorator
39import numpy as np
40from ecmwfapi.api import APIKeyFetchError
41from ecmwfapi.api import APIException as ECMWFAPIException
43from geoarray import GeoArray
44from py_tools_ds.geo.map_info import mapinfo2geotransform
45from pyrsr import RSR
47from ..options.config import GMS_config as CFG
48from . import geoprocessing as GEOP
49from .L1B_P import L1B_object
50from ..model.metadata import get_LayerBandsAssignment
51from ..misc.definition_dicts import get_outFillZeroSaturated, proc_chain, get_mask_classdefinition
52from ..misc.locks import MultiSlotLock
53# from .cloud_masking import Cloud_Mask_Creator # circular dependencies
55from sicor.sicor_ac import ac_gms
56from sicor.sensors import RSImage
57from sicor.Mask import S2Mask
58from sicor.ECMWF import download_variables
60__author__ = 'Daniel Scheffler'
63class L1C_object(L1B_object):
64 def __init__(self, L1B_obj=None):
65 super(L1C_object, self).__init__()
67 if L1B_obj:
68 # populate attributes
69 [setattr(self, key, value) for key, value in L1B_obj.__dict__.items()]
71 # private attributes
72 self._VZA_arr = None
73 self._VAA_arr = None
74 self._SZA_arr = None
75 self._SAA_arr = None
76 self._RAA_arr = None
77 self._lonlat_arr = None
79 self.proc_level = 'L1C'
80 self.proc_status = 'initialized'
82 @property
83 def lonlat_arr(self):
84 """Calculates pixelwise 2D-array with longitude and latitude coordinates.
86 :return:
87 """
88 if self._lonlat_arr is None:
89 self.logger.info('Calculating LonLat array...')
90 self._lonlat_arr = \
91 GEOP.get_lonlat_coord_array(self.shape_fullArr, self.arr_pos,
92 mapinfo2geotransform(self.MetaObj.map_info),
93 self.MetaObj.projection,
94 meshwidth=10, # for faster processing
95 nodata_mask=None, # dont overwrite areas outside the image with nodata
96 outFill=get_outFillZeroSaturated(np.float32)[0])[0]
97 return self._lonlat_arr
99 @lonlat_arr.setter
100 def lonlat_arr(self, lonlat_arr):
101 self._lonlat_arr = lonlat_arr
103 @lonlat_arr.deleter
104 def lonlat_arr(self):
105 self._lonlat_arr = None
107 @property
108 def VZA_arr(self):
109 """Get viewing zenith angle.
111 :return:
112 """
113 if self._VZA_arr is None:
114 self.logger.info('Calculating viewing zenith array...')
115 if self.MetaObj.ViewingAngle_arrProv:
116 # Sentinel-2
117 self._VZA_arr = GEOP.adjust_acquisArrProv_to_shapeFullArr(
118 {k: v.tolist() for k, v in self.MetaObj.ViewingAngle_arrProv.items()},
119 self.shape_fullArr,
120 meshwidth=10, # for faster processing
121 subset=None,
122 bandwise=0)
123 else:
124 self._VZA_arr = GEOP.calc_VZA_array(self.shape_fullArr, self.arr_pos, self.fullSceneCornerPos,
125 float(self.MetaObj.ViewingAngle),
126 float(self.MetaObj.FOV),
127 self.logger,
128 nodata_mask=None, # dont overwrite areas outside image with nodata
129 outFill=get_outFillZeroSaturated(np.float32)[0],
130 meshwidth=10) # for faster processing
131 return self._VZA_arr
133 @VZA_arr.setter
134 def VZA_arr(self, VZA_arr):
135 self._VZA_arr = VZA_arr
137 @VZA_arr.deleter
138 def VZA_arr(self):
139 self._VZA_arr = None
141 @property
142 def VAA_arr(self):
143 """Get viewing azimuth angle.
145 :return:
146 """
147 if self._VAA_arr is None:
148 self.logger.info('Calculating viewing azimuth array...')
149 if self.MetaObj.IncidenceAngle_arrProv:
150 # Sentinel-2
151 self._VAA_arr = GEOP.adjust_acquisArrProv_to_shapeFullArr(
152 {k: v.tolist() for k, v in self.MetaObj.IncidenceAngle_arrProv.items()},
153 self.shape_fullArr,
154 meshwidth=10, # for faster processing
155 subset=None,
156 bandwise=0)
157 else:
158 # only a mean VAA is available
159 if self.VAA_mean is None:
160 self.VAA_mean = \
161 GEOP.calc_VAA_using_fullSceneCornerLonLat(self.fullSceneCornerLonLat, self.MetaObj.orbitParams)
162 assert isinstance(self.VAA_mean, float)
164 self._VAA_arr = np.full(self.VZA_arr.shape, self.VAA_mean, np.float32)
165 return self._VAA_arr
167 @VAA_arr.setter
168 def VAA_arr(self, VAA_arr):
169 self._VAA_arr = VAA_arr
171 @VAA_arr.deleter
172 def VAA_arr(self):
173 self._VAA_arr = None
175 @property
176 def SZA_arr(self):
177 """Get solar zenith angle.
179 :return:
180 """
181 if self._SZA_arr is None:
182 self.logger.info('Calculating solar zenith and azimuth arrays...')
183 self._SZA_arr, self._SAA_arr = \
184 GEOP.calc_SZA_SAA_array(
185 self.shape_fullArr, self.arr_pos,
186 self.MetaObj.AcqDate,
187 self.MetaObj.AcqTime,
188 self.fullSceneCornerPos,
189 self.fullSceneCornerLonLat,
190 self.MetaObj.overpassDurationSec,
191 self.logger,
192 meshwidth=10,
193 nodata_mask=None, # dont overwrite areas outside the image with nodata
194 outFill=get_outFillZeroSaturated(np.float32)[0],
195 accurracy=CFG.SZA_SAA_calculation_accurracy,
196 lonlat_arr=self.lonlat_arr if CFG.SZA_SAA_calculation_accurracy == 'fine' else None)
197 return self._SZA_arr
199 @SZA_arr.setter
200 def SZA_arr(self, SZA_arr):
201 self._SZA_arr = SZA_arr
203 @SZA_arr.deleter
204 def SZA_arr(self):
205 self._SZA_arr = None
207 @property
208 def SAA_arr(self):
209 """Get solar azimuth angle.
211 :return:
212 """
213 if self._SAA_arr is None:
214 # noinspection PyStatementEffect
215 self.SZA_arr # getter also sets self._SAA_arr
216 return self._SAA_arr
218 @SAA_arr.setter
219 def SAA_arr(self, SAA_arr):
220 self._SAA_arr = SAA_arr
222 @SAA_arr.deleter
223 def SAA_arr(self):
224 self._SAA_arr = None
226 @property
227 def RAA_arr(self):
228 """Get relative azimuth angle.
230 :return:
231 """
232 if self._RAA_arr is None:
233 self.logger.info('Calculating relative azimuth array...')
234 self._RAA_arr = GEOP.calc_RAA_array(self.SAA_arr, self.VAA_mean,
235 nodata_mask=None, outFill=get_outFillZeroSaturated(np.float32)[0])
236 return self._RAA_arr
238 @RAA_arr.setter
239 def RAA_arr(self, RAA_arr):
240 self._RAA_arr = RAA_arr
242 @RAA_arr.deleter
243 def RAA_arr(self):
244 self._RAA_arr = None
246 def delete_ac_input_arrays(self):
247 """Delete AC input arrays if they are not needed anymore."""
248 self.logger.info('Deleting input arrays for atmospheric correction...')
249 del self.VZA_arr
250 del self.SZA_arr
251 del self.SAA_arr
252 del self.RAA_arr
253 del self.lonlat_arr
255 # use self.dem deleter
256 # would have to be resampled when writing MGRS tiles
257 # -> better to directly warp it to the output dims and projection
258 del self.dem
261class AtmCorr(object):
262 def __init__(self, *L1C_objs, reporting=False):
263 """Wrapper around atmospheric correction by Andre Hollstein, GFZ Potsdam
265 Creates the input arguments for atmospheric correction from one or multiple L1C_object instance(s) belonging to
266 the same scene ID, performs the atmospheric correction and returns the atmospherically corrected L1C object(s).
268 :param L1C_objs: one or more instances of L1C_object belonging to the same scene ID
269 """
270 # FIXME not yet usable for data < 2012 due to missing ECMWF archive
271 L1C_objs = L1C_objs if isinstance(L1C_objs, tuple) else (L1C_objs,)
273 # hidden attributes
274 self._logger = None
275 self._GSDs = []
276 self._data = {}
277 self._metadata = {}
278 self._nodata = {}
279 self._band_spatial_sampling = {}
280 self._options = {}
282 # assertions
283 scene_IDs = [obj.scene_ID for obj in L1C_objs]
284 assert len(list(set(scene_IDs))) == 1, \
285 "Input GMS objects for 'AtmCorr' must all belong to the same scene ID!. Received %s." % scene_IDs
287 self.inObjs = L1C_objs # type: List[L1C_object]
288 self.reporting = reporting
289 self.ac_input = {} # set by self.run_atmospheric_correction()
290 self.results = None # direct output of external atmCorr module (set by run_atmospheric_correction)
291 self.proc_info = {}
292 self.outObjs = [] # atmospherically corrected L1C objects
294 # append AtmCorr object to input L1C objects
295 # [setattr(L1C_obj, 'AtmCorr', self) for L1C_obj in self.inObjs] # too big for serialization
297 if not re.search(r'Sentinel-2', self.inObjs[0].satellite, re.I):
298 self.logger.debug('Calculation of acquisition geometry arrays is currently only validated for Sentinel-2!')
299 # validation possible by comparing S2 angles provided by ESA with own angles # TODO
301 @property
302 def logger(self):
303 if self._logger and self._logger.handlers[:]:
304 return self._logger
305 else:
306 if len(self.inObjs) == 1:
307 # just use the logger of the inObj
308 logger_atmCorr = self.inObjs[0].logger
309 else:
310 # in case of multiple GMS objects to be processed at once:
311 # get the logger of the first inObj
312 logger_atmCorr = self.inObjs[0].logger
314 # add additional file handlers for the remaining inObj (that belong to the same scene_ID)
315 for inObj in self.inObjs[1:]:
316 path_logfile = inObj.pathGen.get_path_logfile()
317 fileHandler = logging.FileHandler(path_logfile, mode='a')
318 fileHandler.setFormatter(logger_atmCorr.formatter_fileH)
319 fileHandler.setLevel(CFG.log_level)
321 logger_atmCorr.addHandler(fileHandler)
323 inObj.close_loggers()
325 self._logger = logger_atmCorr
326 return self._logger
328 @logger.setter
329 def logger(self, logger):
330 assert isinstance(logger, logging.Logger) or logger in ['not set', None], \
331 "AtmCorr.logger can not be set to %s." % logger
332 if logger in ['not set', None]:
333 self._logger.close()
334 self._logger = logger
335 else:
336 self._logger = logger
338 @logger.deleter
339 def logger(self):
340 if self._logger not in [None, 'not set']:
341 self._logger.close()
342 self._logger = None
344 [inObj.close_loggers() for inObj in self.inObjs]
346 @property
347 def GSDs(self):
348 """
349 Returns a list of spatial samplings within the input GMS objects, e.g. [10,20,60].
350 """
351 for obj in self.inObjs:
352 if obj.arr.xgsd != obj.arr.ygsd:
353 warnings.warn("X/Y GSD is not equal for entity ID %s" % obj.entity_ID +
354 (' (%s)' % obj.subsystem if obj.subsystem else '') +
355 'Using X-GSD as key for spatial sampling dictionary.')
356 self._GSDs.append(obj.arr.xgsd)
358 return self._GSDs
360 @property
361 def data(self):
362 """
364 :return:
365 ___ attribute: data, type:<class 'dict'>
366 ______ key:B05, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 085998540.0803833 ]]
367 ______ key:B01, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 131225590.13208008]]
368 ______ key:B06, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] .14965820.13977051]]
369 ______ key:B11, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] .11492920.10192871]]
370 ______ key:B02, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 104187010.10308838]]
371 ______ key:B10, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 013099670.01300049]]
372 ______ key:B08, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] .16857910.15783691]]
373 ______ key:B04, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 065490720.06228638]]
374 ______ key:B03, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 082702640.08148193]]
375 ______ key:B12, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 068420410.06060791]]
376 ______ key:B8A, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 192138670.17553711]]
377 ______ key:B09, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] .09600830.09887695]]
378 ______ key:B07, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 173339840.15600586]]
379 """
380 if not self._data:
381 data_dict = {}
383 for inObj in self.inObjs:
384 for bandN, bandIdx in inObj.arr.bandnames.items():
385 if bandN not in data_dict:
386 # float32! -> conversion to np.float16 will convert -9999 to -10000
387 arr2pass = inObj.arr[:, :, bandIdx].astype(np.float32)
388 arr2pass[arr2pass == inObj.arr.nodata] = np.nan # set nodata values to np.nan
389 data_dict[bandN] = (arr2pass / inObj.MetaObj.ScaleFactor).astype(np.float16)
390 else:
391 inObj.logger.warning("Band '%s' cannot be included into atmospheric correction because it "
392 "exists multiple times." % bandN)
394 # validate: data must have all bands needed for AC
395 full_LBA = get_LayerBandsAssignment(self.inObjs[0].GMS_identifier, return_fullLBA=True)
396 all_bNs_AC = ['B%s' % i if len(i) == 2 else 'B0%s' % i for i in full_LBA]
397 if not all([bN in list(data_dict.keys()) for bN in all_bNs_AC]):
398 raise RuntimeError('Atmospheric correction did not receive all the needed bands. \n\tExpected: %s;\n\t'
399 'Received: %s' % (str(all_bNs_AC), str(list(sorted(data_dict.keys())))))
401 self._data = data_dict
403 return self._data
405 @data.setter
406 def data(self, data_dict):
407 assert isinstance(data_dict, dict), \
408 "'data' can only be set to a dictionary with band names as keys and numpy arrays as values."
409 self._data = data_dict
411 @property
412 def nodata(self):
413 """
415 :return:
416 ___ attribute: nodata, type:<class 'dict'>
417 ______ key:60.0, value_type:<class 'numpy.ndarray'>, repr: [[ TrueTrueTrue ..., [..] False False False]]
418 ______ key:10.0, value_type:<class 'numpy.ndarray'>, repr: [[ TrueTrueTrue ..., [..] False False False]]
419 ______ key:20.0, value_type:<class 'numpy.ndarray'>, repr: [[ TrueTrueTrue ..., [..] False False False]]
420 """
422 if not self._nodata:
423 for inObj in self.inObjs:
424 self._nodata[inObj.arr.xgsd] = ~inObj.arr.mask_nodata[:]
426 return self._nodata
428 @property
429 def tile_name(self):
430 """Returns S2A tile name.
431 NOTE: this is only needed if no DEM is passed to ac_gms
433 :return: e.g.
434 '32UMA'
435 """
437 return '' # FIXME
439 @property
440 def band_spatial_sampling(self):
441 """
443 :return: e.g.
444 {'B01': 60.0,
445 'B02': 10.0,
446 'B03': 10.0,
447 'B04': 10.0,
448 'B05': 20.0,
449 'B06': 20.0,
450 'B07': 20.0,
451 'B08': 10.0,
452 'B09': 60.0,
453 'B10': 60.0,
454 'B11': 20.0,
455 'B12': 20.0,
456 'B8A': 20.0}
457 """
459 if not self._band_spatial_sampling:
460 for inObj in self.inObjs:
461 for bandN in inObj.arr.bandnames:
462 if bandN not in self._band_spatial_sampling:
463 self._band_spatial_sampling[bandN] = inObj.arr.xgsd
464 return self._band_spatial_sampling
466 @property
467 def metadata(self):
468 """
470 :return:
471 ___ attribute: metadata, type:<class 'dict'>
472 ______ key:spatial_samplings
473 _________ key:60.0
474 ____________ key:ULY, value_type:<class 'int'>, repr: 4900020
475 ____________ key:NCOLS, value_type:<class 'int'>, repr: 1830
476 ____________ key:XDIM, value_type:<class 'int'>, repr: 60
477 ____________ key:ULX, value_type:<class 'int'>, repr: 600000
478 ____________ key:NROWS, value_type:<class 'int'>, repr: 1830
479 ____________ key:YDIM, value_type:<class 'int'>, repr: -60
480 _________ key:10.0
481 ____________ key:ULY, value_type:<class 'int'>, repr: 4900020
482 ____________ key:NCOLS, value_type:<class 'int'>, repr: 10980
483 ____________ key:XDIM, value_type:<class 'int'>, repr: 10
484 ____________ key:ULX, value_type:<class 'int'>, repr: 600000
485 ____________ key:NROWS, value_type:<class 'int'>, repr: 10980
486 ____________ key:YDIM, value_type:<class 'int'>, repr: -10
487 _________ key:20.0
488 ____________ key:ULY, value_type:<class 'int'>, repr: 4900020
489 ____________ key:NCOLS, value_type:<class 'int'>, repr: 5490
490 ____________ key:XDIM, value_type:<class 'int'>, repr: 20
491 ____________ key:ULX, value_type:<class 'int'>, repr: 600000
492 ____________ key:NROWS, value_type:<class 'int'>, repr: 5490
493 ____________ key:YDIM, value_type:<class 'int'>, repr: -20
494 ______ key:SENSING_TIME, value_type:<class 'datetime.datetime'>, repr: 2016-03-26 10:34:06.538000+00:00
495 """
497 if not self._metadata:
498 del self.logger # otherwise each input object would have multiple fileHandlers
500 metadata = dict(
501 U=self.inObjs[0].MetaObj.EarthSunDist,
502 SENSING_TIME=self.inObjs[0].acq_datetime,
503 # SENSING_TIME=datetime.strptime('2015-08-12 10:40:21 +0000', '%Y-%m-%d %H:%M:%S %z'),
504 viewing_zenith=self._meta_get_viewing_zenith(),
505 viewing_azimuth=self._meta_get_viewing_azimuth(),
506 relative_viewing_azimuth=self._meta_get_relative_viewing_azimuth(),
507 sun_mean_azimuth=self.inObjs[0].MetaObj.SunAzimuth,
508 sun_mean_zenith=90 - self.inObjs[0].MetaObj.SunElevation,
509 solar_irradiance=self._meta_get_solar_irradiance(),
510 aux_data=self._meta_get_aux_data(),
511 spatial_samplings=self._meta_get_spatial_samplings()
512 )
514 self._metadata = metadata
516 return self._metadata
518 @property
519 def options(self):
520 # type: () -> dict
521 """Returns a dictionary containing AC options.
522 """
523 if self._options:
524 return self._options
525 else:
526 self._options = self.inObjs[0].ac_options
527 self._options["AC"]['bands'] = [b for b in self.data.keys() if b in self._options["AC"]['bands']]
528 self._options["report"]["reporting"] = self.reporting
529 return self._options
531 def _meta_get_spatial_samplings(self):
532 """
534 :return:
535 {10.0: {'NCOLS': 10980,
536 'NROWS': 10980,
537 'ULX': 499980.0,
538 'ULY': 5800020.0,
539 'XDIM': 10.0,
540 'YDIM': -10.0},
541 20.0: {'NCOLS': 5490,
542 'NROWS': 5490,
543 'ULX': 499980.0,
544 'ULY': 5800020.0,
545 'XDIM': 20.0,
546 'YDIM': -20.0},
547 60.0: {'NCOLS': 1830,
548 'NROWS': 1830,
549 'ULX': 499980.0,
550 'ULY': 5800020.0,
551 'XDIM': 60.0,
552 'YDIM': -60.0}}
553 """
554 # set corner coordinates and dims
555 spatial_samplings = {}
557 for inObj in self.inObjs:
559 # validate GSD
560 if inObj.arr.xgsd != inObj.arr.ygsd:
561 warnings.warn("X/Y GSD is not equal for entity ID %s" % inObj.entity_ID +
562 (' (%s)' % inObj.subsystem if inObj.subsystem else '') +
563 'Using X-GSD as key for spatial sampling dictionary.')
565 # set spatial information
566 spatial_samplings[inObj.arr.xgsd] = dict(
567 ULX=inObj.arr.box.boxMapYX[0][1],
568 ULY=inObj.arr.box.boxMapYX[0][0],
569 XDIM=inObj.arr.xgsd,
570 YDIM=-inObj.arr.ygsd,
571 NROWS=inObj.arr.rows,
572 NCOLS=inObj.arr.cols)
574 return spatial_samplings
576 def _meta_get_solar_irradiance(self):
577 """
579 :return:
580 {'B01': 1913.57,
581 'B02': 1941.63,
582 'B03': 1822.61,
583 'B04': 1512.79,
584 'B05': 1425.56,
585 'B06': 1288.32,
586 'B07': 1163.19,
587 'B08': 1036.39,
588 'B09': 813.04,
589 'B10': 367.15,
590 'B11': 245.59,
591 'B12': 85.25,
592 'B8A': 955.19}
593 """
595 solar_irradiance = {}
597 for inObj in self.inObjs:
598 for bandN in inObj.arr.bandnames:
599 lba_key = bandN[2:] if bandN.startswith('B0') else bandN[1:]
600 if bandN not in solar_irradiance:
601 solar_irradiance[bandN] = inObj.MetaObj.SolIrradiance[lba_key]
603 return solar_irradiance
605 def _meta_get_viewing_zenith(self):
606 """
608 :return: {B10:ndarray(dtype=float16),[...],B09:ndarray(dtype=float16)}
609 """
611 viewing_zenith = {}
613 for inObj in self.inObjs: # type: L1C_object
614 for bandN, bandIdx in inObj.arr.bandnames.items():
615 if bandN not in viewing_zenith:
616 arr2pass = inObj.VZA_arr[:, :, bandIdx] if inObj.VZA_arr.ndim == 3 else inObj.VZA_arr
617 viewing_zenith[bandN] = arr2pass.astype(np.float16)
618 # viewing_zenith[bandN] = inObj.VZA_arr[:, :, bandIdx] if inObj.VZA_arr.ndim==3 else inObj.VZA_arr
619 return viewing_zenith
621 def _meta_get_viewing_azimuth(self):
622 """
624 :return: {B10:ndarray(dtype=float16),[...],B09:ndarray(dtype=float16)}
625 """
627 viewing_azimuth = {}
629 for inObj in self.inObjs: # type: L1C_object
630 for bandN, bandIdx in inObj.arr.bandnames.items():
631 if bandN not in viewing_azimuth:
632 arr2pass = inObj.VAA_arr[:, :, bandIdx] if inObj.VAA_arr.ndim == 3 else inObj.VAA_arr
633 viewing_azimuth[bandN] = arr2pass.astype(np.float16)
634 # viewing_azimuth[bandN] = inObj.VAA_arr[:, :, bandIdx] if inObj.VAA_arr.ndim==3 else inObj.VAA_arr
636 return viewing_azimuth
638 def _meta_get_relative_viewing_azimuth(self):
639 """
641 :return: {B10:ndarray(dtype=float16),[...],B09:ndarray(dtype=float16)}
642 """
644 relative_viewing_azimuth = {}
646 for inObj in self.inObjs: # type: L1C_object
647 for bandN, bandIdx in inObj.arr.bandnames.items():
648 if bandN not in relative_viewing_azimuth:
649 arr2pass = inObj.RAA_arr[:, :, bandIdx] if inObj.RAA_arr.ndim == 3 else inObj.RAA_arr
650 relative_viewing_azimuth[bandN] = arr2pass.astype(np.float16)
651 # relative_viewing_azimuth[bandN] = \
652 # inObj.RAA_arr[:, :, bandIdx] if inObj.RAA_arr.ndim==3 else inObj.RAA_arr
654 return relative_viewing_azimuth
656 def _meta_get_aux_data(self):
657 """
659 :return: {lons:ndarray(dtype=float16),,lats:ndarray(dtype=float16)}
660 """
662 aux_data = dict(
663 # set lons and lats (a 2D array for all bands is enough (different band resolutions dont matter))
664 lons=self.inObjs[0].lonlat_arr[::10, ::10, 0].astype(np.float16), # 2D array of lon values: 0° - 360°
665 lats=self.inObjs[0].lonlat_arr[::10, ::10, 1].astype(np.float16) # 2D array of lat values: -90° - 90°
666 # FIXME correct to reduce resolution here by factor 10?
667 )
669 return aux_data
671 def _get_dem(self):
672 """Get a DEM to be used in atmospheric correction.
674 :return: <np.ndarray> 2D array (with 20m resolution in case of Sentinel-2)
675 """
676 # determine which input GMS object is used to generate DEM
677 if re.search(r'Sentinel-2', self.inObjs[0].satellite):
678 # in case of Sentinel-2 the 20m DEM must be passed
679 inObj4dem = [obj for obj in self.inObjs if obj.arr.xgsd == 20]
680 if not inObj4dem:
681 self.logger.warning('Sentinel-2 20m subsystem could not be found. DEM passed to '
682 'atmospheric correction might have wrong resolution.')
683 inObj4dem = inObj4dem[0]
684 else:
685 inObj4dem = self.inObjs[0]
687 try:
688 dem = inObj4dem.dem[:].astype(np.float32)
689 except Exception as e:
690 dem = None
691 self.logger.warning('A static elevation is assumed during atmospheric correction due to an error during '
692 'creation of the DEM corresponding to scene %s (entity ID: %s). Error message was: '
693 '\n%s\n' % (self.inObjs[0].scene_ID, self.inObjs[0].entity_ID, repr(e)))
694 self.logger.info("Print traceback in case you care:")
695 self.logger.warning(traceback.format_exc())
697 return dem
699 def _get_srf(self):
700 """Returns an instance of SRF in the same structure like sicor.sensors.SRF.SensorSRF
701 """
702 # FIXME calculation of center wavelengths within RSR() used not the GMS algorithm
703 # SRF instance must be created for all bands and the previous proc level
704 GMSid_fullScene = self.inObjs[0].GMS_identifier
705 GMSid_fullScene.subsystem = ''
706 GMSid_fullScene.proc_level = proc_chain[proc_chain.index(self.inObjs[0].proc_level) - 1]
708 return RSR(satellite=GMSid_fullScene.satellite,
709 sensor=GMSid_fullScene.sensor,
710 subsystem=GMSid_fullScene.subsystem,
711 wvl_unit='nanometers',
712 format_bandnames=True,
713 no_pan=CFG.skip_pan,
714 no_thermal=CFG.skip_thermal,
715 after_ac=GMSid_fullScene.proc_level in ['L1C', 'L2A', 'L2B', 'L2C'],
716 sort_by_cwl=CFG.sort_bands_by_cwl
717 )
719 def _get_mask_clouds(self):
720 """Returns an instance of S2Mask in case cloud mask is given by input GMS objects. Otherwise None is returned.
722 :return:
723 """
725 tgt_res = self.inObjs[0].ac_options['cld_mask']['target_resolution']
727 # check if input GMS objects provide a cloud mask
728 avail_cloud_masks = {inObj.GMS_identifier.subsystem: inObj.mask_clouds for inObj in self.inObjs}
729 no_avail_CMs = list(set(avail_cloud_masks.values())) == [None]
731 # compute cloud mask if not already provided
732 if no_avail_CMs:
733 algorithm = CFG.cloud_masking_algorithm[self.inObjs[0].satellite]
735 if algorithm == 'SICOR':
736 return None
738 else:
739 # FMASK or Classical Bayesian
740 try:
741 from .cloud_masking import Cloud_Mask_Creator
743 CMC = Cloud_Mask_Creator(self.inObjs[0], algorithm=algorithm, tempdir_root=CFG.path_tempdir)
744 CMC.calc_cloud_mask()
745 cm_geoarray = CMC.cloud_mask_geoarray
746 cm_array = CMC.cloud_mask_array
747 cm_legend = CMC.cloud_mask_legend
748 except Exception:
749 self.logger.error('\nAn error occurred during FMASK cloud masking. Error message was: ')
750 self.logger.error(traceback.format_exc())
751 return None
753 else:
754 # check if there is a cloud mask with suitable GSD
755 inObjs2use = [obj for obj in self.inObjs if obj.mask_clouds is not None and obj.mask_clouds.xgsd == tgt_res]
756 if not inObjs2use:
757 raise ValueError('Error appending cloud mask to input arguments of atmospheric correction. No input '
758 'GMS object provides a cloud mask with spatial resolution of %s.' % tgt_res)
759 inObj2use = inObjs2use[0]
761 # get mask (geo)array
762 cm_geoarray = inObj2use.mask_clouds
763 cm_array = inObj2use.mask_clouds[:]
765 # get legend
766 cm_legend = get_mask_classdefinition('mask_clouds', inObj2use.satellite)
767 # {'Clear': 10, 'Thick Clouds': 20, 'Thin Clouds': 30, 'Snow': 40} # FIXME hardcoded
769 # validate that xGSD equals yGSD
770 if cm_geoarray.xgsd != cm_geoarray.ygsd:
771 warnings.warn("Cloud mask X/Y GSD is not equal for entity ID %s" % inObj2use.entity_ID +
772 (' (%s)' % inObj2use.subsystem if inObj2use.subsystem else '') +
773 'Using X-GSD as key for cloud mask geocoding.')
775 # get geocoding
776 cm_geocoding = self.metadata["spatial_samplings"][tgt_res]
778 # get nodata value
779 self.options['cld_mask']['nodata_value_mask'] = cm_geoarray.nodata
781 # append cloud mask to input object with the same spatial resolution if there was no mask before
782 for inObj in self.inObjs:
783 if inObj.arr.xgsd == cm_geoarray.xgsd:
784 inObj.mask_clouds = cm_geoarray
785 inObj.build_combined_masks_array()
786 break # appending it to one inObj is enough
788 return S2Mask(mask_array=cm_array,
789 mask_legend=cm_legend,
790 geo_coding=cm_geocoding)
792 def _check_or_download_ECMWF_data(self):
793 """Check if ECMWF files are already downloaded. If not, start the downloader."""
794 self.logger.info('Checking if ECMWF data are available... (if not, run download!)')
796 default_products = [
797 "fc_T2M",
798 "fc_O3",
799 "fc_SLP",
800 "fc_TCWV",
801 "fc_GMES_ozone",
802 "fc_total_AOT_550nm",
803 "fc_sulphate_AOT_550nm",
804 "fc_black_carbon_AOT_550nm",
805 "fc_dust_AOT_550nm",
806 "fc_organic_matter_AOT_550nm",
807 "fc_sea_salt_AOT_550nm"]
809 # NOTE: use_signals must be set to True while executed as multiprocessing worker (e.g., in multiprocessing.Pool)
810 @timeout_decorator.timeout(seconds=60*10, timeout_exception=TimeoutError)
811 def run_request():
812 try:
813 with MultiSlotLock('ECMWF download lock', allowed_slots=1, logger=self.logger):
814 t0 = time()
815 # NOTE: download_variables does not accept a logger -> so the output may be invisible in WebApp
816 results = download_variables(date_from=self.inObjs[0].acq_datetime,
817 date_to=self.inObjs[0].acq_datetime,
818 db_path=CFG.path_ECMWF_db,
819 max_step=120, # default
820 ecmwf_variables=default_products,
821 processes=0, # singleprocessing
822 force=False) # dont force download if files already exist
823 t1 = time()
824 self.logger.info("Runtime: %.2f" % (t1 - t0))
825 for result in results:
826 self.logger.info(result)
828 except APIKeyFetchError:
829 self.logger.error("ECMWF data download failed due to missing API credentials.")
831 except (ECMWFAPIException, Exception):
832 self.logger.error("ECMWF data download failed for scene %s (entity ID: %s). Traceback: "
833 % (self.inObjs[0].scene_ID, self.inObjs[0].entity_ID))
834 self.logger.error(traceback.format_exc())
836 try:
837 run_request()
838 except TimeoutError:
839 self.logger.error("ECMWF data download failed due to API request timeout after waiting 5 minutes.")
841 def _validate_snr_source(self):
842 """Check if the given file path for the SNR model exists - if not, use a constant SNR of 500."""
843 if not os.path.isfile(self.options["uncertainties"]["snr_model"]):
844 self.logger.warning('No valid SNR model found for %s %s. Using constant SNR to compute uncertainties of '
845 'atmospheric correction.' % (self.inObjs[0].satellite, self.inObjs[0].sensor))
846 # self.options["uncertainties"]["snr_model"] = np.nan # causes the computed uncertainties to be np.nan
847 self.options["uncertainties"]["snr_model"] = 500 # use a constant SNR of 500 to compute uncertainties
849 def run_atmospheric_correction(self, dump_ac_input=False):
850 # type: (bool) -> list
851 """Collects all input data for atmospheric correction, runs the AC and returns the corrected L1C objects
852 containing surface reflectance.
854 :param dump_ac_input: allows to dump the inputs of AC to the scene's processing folder in case AC fails
855 :return: list of L1C_object instances containing atmospherically corrected data
856 """
858 # collect input args/kwargs for AC
859 self.logger.info('Calculating input data for atmospheric correction...')
861 rs_data = dict(
862 data=self.data,
863 metadata=self.metadata,
864 nodata=self.nodata,
865 band_spatial_sampling=self.band_spatial_sampling,
866 tile_name=self.tile_name,
867 dem=self._get_dem(),
868 srf=self._get_srf(),
869 mask_clouds=self._get_mask_clouds()
870 # returns an instance of S2Mask or None if cloud mask is not given by input GMS objects
871 ) # NOTE: all keys of this dict are later converted to attributes of RSImage
873 # remove empty values from RSImage kwargs because SICOR treats any kind of RSImage attributes as given
874 # => 'None'-attributes may cause issues
875 rs_data = {k: v for k, v in rs_data.items() if v is not None}
877 script = False
879 # check if ECMWF data are available - if not, start the download
880 if CFG.auto_download_ecmwf:
881 self._check_or_download_ECMWF_data()
883 # validate SNR
884 if CFG.ac_estimate_accuracy:
885 self._validate_snr_source()
887 # create an instance of RSImage
888 rs_image = RSImage(**rs_data)
890 self.ac_input = dict(
891 rs_image=rs_image,
892 options=self.options, # dict
893 logger=repr(self.logger), # only a string
894 script=script
895 )
897 # path_dump = self.inObjs[0].pathGen.get_path_ac_input_dump()
898 # with open(path_dump, 'wb') as outF:
899 # dill.dump(self.ac_input, outF)
901 # run AC
902 self.logger.info('Atmospheric correction started.')
903 try:
904 rs_image.logger = self.logger
905 self.results = ac_gms(rs_image, self.options, logger=self.logger, script=script)
907 except Exception as e:
908 self.logger.error('\nAn error occurred during atmospheric correction. BE AWARE THAT THE SCENE %s '
909 '(ENTITY ID %s) HAS NOT BEEN ATMOSPHERICALLY CORRECTED! Error message was: \n%s\n'
910 % (self.inObjs[0].scene_ID, self.inObjs[0].entity_ID, repr(e)))
911 self.logger.error(traceback.format_exc())
912 # TODO include that in the job summary
914 # serialize AC input
915 if dump_ac_input:
916 path_dump = self.inObjs[0].pathGen.get_path_ac_input_dump()
917 with open(path_dump, 'wb') as outF:
918 self.ac_input['rs_image'].logger = None
919 dill.dump(self.ac_input, outF)
921 self.logger.error('An error occurred during atmospheric correction. Inputs have been dumped to %s.'
922 % path_dump)
924 # delete AC input arrays
925 for inObj in self.inObjs: # type: L1C_object
926 inObj.delete_ac_input_arrays()
928 return list(self.inObjs)
930 finally:
931 # rs_image.logger must be closed properly in any case
932 if rs_image.logger is not None:
933 rs_image.logger.close()
935 # get processing infos
936 self.proc_info = self.ac_input['options']['processing']
938 # join results
939 self._join_results_to_inObjs() # sets self.outObjs
941 # delete input arrays that are not needed anymore
942 [inObj.delete_ac_input_arrays() for inObj in self.inObjs]
944 return self.outObjs
946 def _join_results_to_inObjs(self):
947 """
948 Join results of atmospheric correction to the input GMS objects.
949 """
951 self.logger.info('Joining results of atmospheric correction to input GMS objects.')
952 # delete logger
953 # -> otherwise logging in inObjs would open a second FileHandler to the same file (which is permitted)
954 del self.logger
956 self._join_data_ac()
957 self._join_data_errors(bandwise=CFG.ac_bandwise_accuracy)
958 self._join_mask_clouds()
959 self._join_mask_confidence_array()
961 # update masks (always do that because masks can also only contain one layer)
962 [inObj.build_combined_masks_array() for inObj in self.inObjs]
964 # update AC processing info
965 [inObj.ac_options['processing'].update(self.proc_info) for inObj in self.inObjs]
967 self.outObjs = self.inObjs
969 def _join_data_ac(self):
970 """
971 Join ATMOSPHERICALLY CORRECTED ARRAY as 3D int8 or int16 BOA reflectance array, scaled to scale factor from
972 config.
973 """
975 if self.results.data_ac is not None:
976 for inObj in self.inObjs:
977 self.logger.info('Joining image data to %s.' % (inObj.entity_ID if not inObj.subsystem else
978 '%s %s' % (inObj.entity_ID, inObj.subsystem)))
980 assert isinstance(inObj, L1B_object)
981 nodata = self.results.nodata[inObj.arr.xgsd] # 2D mask with True outside of image coverage
982 ac_bandNs = [bandN for bandN in inObj.arr.bandnames if bandN in self.results.data_ac.keys()]
983 out_LBA = [bN.split('B0')[1] if bN.startswith('B0') else bN.split('B')[1] for bN in ac_bandNs]
985 # update metadata #
986 ###################
988 inObj.arr_desc = 'BOA_Ref'
989 inObj.MetaObj.bands = len(self.results.data_ac)
990 inObj.MetaObj.PhysUnit = 'BOA_Reflectance in [0-%d]' % CFG.scale_factor_BOARef
991 inObj.MetaObj.LayerBandsAssignment = out_LBA
992 inObj.LayerBandsAssignment = out_LBA
993 inObj.MetaObj.filter_layerdependent_metadata()
995 ##################################################################################
996 # join SURFACE REFLECTANCE as 3D int16 array, scaled to scale factor from config #
997 ##################################################################################
999 oF_refl, oZ_refl, oS_refl = get_outFillZeroSaturated(inObj.arr.dtype)
1000 surf_refl = np.dstack([self.results.data_ac[bandN] for bandN in ac_bandNs])
1001 surf_refl *= CFG.scale_factor_BOARef # scale using scale factor (output is float16)
1003 # set AC nodata values to GMS outFill
1004 # NOTE: AC nodata contains a pixel mask where at least one band is no data
1005 # => Setting these pixels to outZero would also reduce pixel values of surrounding pixels in
1006 # spatial homogenization (because resampling only ignores -9999).
1007 # It would be possible to generate a zero-data mask here for each subsystem and apply it after
1008 # spatial homogenization. Alternatively zero-data pixels could be interpolated spectrally or
1009 # spatially within L1A processor (also see issue #74).
1010 surf_refl[nodata] = oF_refl # overwrite AC nodata values with GMS outFill
1012 # apply the original nodata mask (indicating background values)
1013 surf_refl[np.array(inObj.mask_nodata).astype(np.int8) == 0] = oF_refl
1015 # set AC NaNs to GMS outFill
1016 # NOTE: SICOR result has NaNs at no data positions AND non-clear positions
1017 if self.results.bad_data_value is np.nan:
1018 surf_refl[np.isnan(surf_refl)] = oF_refl
1019 else:
1020 surf_refl[surf_refl == self.results.bad_data_value] = oF_refl
1022 # use inObj.arr setter to generate a GeoArray
1023 inObj.arr = surf_refl.astype(inObj.arr.dtype) # -> int16 (also converts NaNs to 0 if needed
1025 else:
1026 self.logger.warning('Atmospheric correction did not return a result for the input array. '
1027 'Thus the output keeps NOT atmospherically corrected.')
1029 def _join_data_errors(self, bandwise=False):
1030 """Join ERRORS ARRAY as 3D or 2D int8 or int16 BOA reflectance array, scaled to scale factor from config.
1032 :param bandwise: if True, a 3D array with bandwise information for each pixel is generated
1033 :return:
1034 """
1035 # TODO ac_error values are very close to 0 -> a scale factor of 255 yields int8 values below 10
1036 # TODO => better to stretch the whole array to values between 0 and 100 and save original min/max?
1037 if self.results.data_errors is not None:
1039 for inObj in self.inObjs:
1040 inObj.logger.info('Joining AC errors to %s.' % (inObj.entity_ID if not inObj.subsystem else
1041 '%s %s' % (inObj.entity_ID, inObj.subsystem)))
1043 nodata = self.results.nodata[inObj.arr.xgsd] # 2D mask with True outside of image coverage
1044 ac_bandNs = [bandN for bandN in inObj.arr.bandnames if bandN in self.results.data_ac.keys()]
1045 out_dtype = np.int8 if CFG.ac_scale_factor_errors <= 255 else np.int16
1046 out_nodata_val = get_outFillZeroSaturated(out_dtype)[0]
1048 # generate raw ac_errors array
1049 ac_errors = np.dstack([self.results.data_errors[bandN] for bandN in ac_bandNs])
1051 # apply scale factor from config to data pixels and overwrite nodata area with nodata value
1052 ac_errors *= CFG.ac_scale_factor_errors # scale using scale factor (output is float16)
1053 ac_errors[np.isnan(ac_errors)] = out_nodata_val # replace NaNs with outnodata
1054 ac_errors[nodata] = out_nodata_val # set all positions where SICOR nodata mask is 1 to outnodata
1055 ac_errors = np.around(ac_errors).astype(out_dtype) # round floats to next even int8/int16 value
1057 # average array over bands if no bandwise information is requested
1058 if not bandwise:
1059 # in case of only one subsystem: directly compute median errors here
1060 if len(self.inObjs) == 1:
1061 ac_errors = np.median(ac_errors, axis=2).astype(ac_errors.dtype)
1063 # in case of multiple subsystems: dont compute median here but first apply geometric homogenization
1064 # -> median could only be computed for each subsystem separately
1065 else:
1066 pass
1068 # set inObj.ac_errors
1069 inObj.ac_errors = ac_errors # setter generates a GeoArray with the same bandnames like inObj.arr
1071 elif not CFG.ac_estimate_accuracy:
1072 self.logger.info("Atmospheric correction did not provide a 'data_errors' array because "
1073 "'ac_estimate_accuracy' was set to False in the job configuration.")
1074 else:
1075 self.logger.warning("Atmospheric correction did not provide a 'data_errors' array. Maybe due to "
1076 "missing SNR model? GMS_object.ac_errors kept None.")
1078 def _join_mask_clouds(self):
1079 """
1080 Join CLOUD MASK as 2D uint8 array.
1081 NOTE: mask_clouds has also methods 'export_mask_rgb()', 'export_confidence_to_jpeg2000()', ...
1082 """
1084 if self.results.mask_clouds.mask_array is not None:
1085 mask_clouds_ac = self.results.mask_clouds.mask_array # uint8 2D array
1087 joined = False
1088 for inObj in self.inObjs:
1090 # delete all previous cloud masks
1091 del inObj.mask_clouds # FIXME validate if FMask product is within AC results
1093 # append mask_clouds only to the input GMS object with the same dimensions
1094 if inObj.arr.shape[:2] == mask_clouds_ac.shape:
1095 inObj.logger.info('Joining mask_clouds to %s.' % (inObj.entity_ID if not inObj.subsystem else
1096 '%s %s' % (inObj.entity_ID, inObj.subsystem)))
1098 inObj.mask_clouds = mask_clouds_ac
1099 inObj.mask_clouds.legend = self.results.mask_clouds.mask_legend # dict(value=string, string=value))
1100 # FIXME legend is not used later
1102 # set cloud mask nodata value
1103 tgt_nodata = get_outFillZeroSaturated(mask_clouds_ac.dtype)[0]
1104 ac_out_nodata = self.ac_input['options']['cld_mask']['nodata_value_mask']
1105 if tgt_nodata not in self.results.mask_clouds.mask_legend.keys():
1106 inObj.mask_clouds[inObj.mask_clouds[:] == ac_out_nodata] = tgt_nodata
1107 mask_clouds_nodata = tgt_nodata
1108 else:
1109 warnings.warn('The cloud mask from AC output already uses the desired nodata value %s for the '
1110 'class %s. Using AC output nodata value %s.'
1111 % (tgt_nodata, self.results.mask_clouds.mask_legend[tgt_nodata], ac_out_nodata))
1112 mask_clouds_nodata = ac_out_nodata
1114 inObj.mask_clouds.nodata = mask_clouds_nodata
1116 joined = True
1118 if not joined:
1119 self.logger.warning('Cloud mask has not been appended to one of the AC inputs because there was no'
1120 'input GMS object with the same dimensions.')
1122 else:
1123 self.logger.warning("Atmospheric correction did not provide a 'mask_clouds.mask_array' array. "
1124 "GMS_object.mask_clouds kept None.")
1126 def _join_mask_confidence_array(self):
1127 """
1128 Join confidence array for mask_clouds.
1129 """
1130 if self.results.mask_clouds.mask_confidence_array is not None and CFG.ac_estimate_accuracy:
1131 cfd_arr = self.results.mask_clouds.mask_confidence_array # float32 2D array, scaled [0-1, nodata 255]
1132 cfd_arr[cfd_arr == self.ac_input['options']['cld_mask']['nodata_value_mask']] = -1
1133 cfd_arr = (cfd_arr * CFG.scale_factor_BOARef).astype(np.int16)
1134 cfd_arr[cfd_arr == -CFG.scale_factor_BOARef] = get_outFillZeroSaturated(cfd_arr.dtype)[0]
1136 joined = False
1137 for inObj in self.inObjs:
1139 # append mask_clouds confidence only to the input GMS object with the same dimensions
1140 if inObj.arr.shape[:2] == cfd_arr.shape:
1141 inObj.logger.info(
1142 'Joining mask_clouds_confidence to %s.' % (inObj.entity_ID if not inObj.subsystem else
1143 '%s %s' % (inObj.entity_ID, inObj.subsystem)))
1145 # set cloud mask confidence array
1146 inObj.mask_clouds_confidence = GeoArray(cfd_arr, inObj.arr.gt, inObj.arr.prj,
1147 nodata=get_outFillZeroSaturated(cfd_arr.dtype)[0])
1148 joined = True
1150 if not joined:
1151 self.logger.warning('Cloud mask confidence array has not been appended to one of the AC inputs because '
1152 'there was no input GMS object with the same dimensions.')
1154 elif CFG.cloud_masking_algorithm[self.inObjs[0].satellite] == 'FMASK':
1155 self.logger.info("Cloud mask confidence array is not appended to AC outputs because "
1156 "'cloud_masking_algorithm' was set to FMASK in the job configuration and FMASK does not "
1157 "confidence arrays.")
1159 elif not CFG.ac_estimate_accuracy:
1160 self.logger.info("Cloud mask confidence array is not appended to AC outputs because "
1161 "'ac_estimate_accuracy' was set to False in the job configuration.")
1162 else:
1163 self.logger.warning("Atmospheric correction did not provide a 'mask_confidence_array' array for "
1164 "attribute 'mask_clouds. GMS_object.mask_clouds_confidence kept None.")