Commit 92ef7f6e authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Added first running version out accuracy layers + writers.

Added options 'ac_bandwise_accuracy', ''spechomo_bandwise_accuracy'.
Moved mask_clouds_confidence, ac_errors, spec_homo_errors to GMS_object.
parent 42865754
......@@ -910,7 +910,7 @@ class AtmCorr(object):
del self.logger
self._join_data_ac()
self._join_data_errors()
self._join_data_errors(bandwise=CFG.ac_bandwise_accuracy)
self._join_mask_clouds()
self._join_mask_confidence_array()
......@@ -967,23 +967,34 @@ class AtmCorr(object):
self.logger.warning('Atmospheric correction did not return a result for the input array. '
'Thus the output keeps NOT atmospherically corrected.')
def _join_data_errors(self):
"""
Join ERRORS ARRAY as 3D int8 or int16 BOA reflectance array, scaled to scale factor from config.
"""
def _join_data_errors(self, bandwise=False):
"""Join ERRORS ARRAY as 3D or 2D int8 or int16 BOA reflectance array, scaled to scale factor from config.
:param bandwise: if True, a 3D array with bandwise information for each pixel is generated
:return:
"""
# TODO how to handle nans?
if self.results.data_errors is not None:
for inObj in self.inObjs:
nodata = self.results.nodata[inObj.arr.xgsd] # 2D mask with True outside of image coverage
ac_bandNs = [bandN for bandN in inObj.arr.bandnames if bandN in self.results.data_ac.keys()]
# generate raw ac_errors array
ac_errors = np.dstack((self.results.data_errors[bandN] for bandN in ac_bandNs))
if not bandwise:
ac_errors = np.nanmedian(ac_errors, axis=2).astype(ac_errors.dtype)
# apply scale factor from config to data pixels and overwrite nodata area with nodata value
ac_errors *= CFG.ac_scale_factor_errors # scale using scale factor (output is float16)
out_dtype = np.int8 if CFG.ac_scale_factor_errors <= 255 else np.int16
ac_errors[nodata] = get_outFillZeroSaturated(out_dtype)[0]
ac_errors = ac_errors.astype(out_dtype)
# set inObj.ac_errors
inObj.ac_errors = ac_errors # setter generates a GeoArray with the same bandnames like inObj.arr
# TODO how to handle nans?
elif not CFG.ac_estimate_accuracy:
self.logger.info("Atmospheric correction did not provide a 'data_errors' array because "
"'ac_estimate_accuracy' was set to False in the job configuration.")
......
......@@ -104,6 +104,7 @@ class L2B_object(L2A_object):
tgt_LBA=tgt_LBA,
nodataVal=self.arr.nodata,
compute_errors=CFG.spechomo_estimate_accuracy,
bandwise_errors=CFG.spechomo_bandwise_accuracy,
fallback_argskwargs=dict(
args=dict(source_CWLs=src_cwls, target_CWLs=tgt_cwls,),
kwargs=dict(kind='linear')
......@@ -162,7 +163,7 @@ class SpectralHomogenizer(object):
def predict_by_machine_learner(self, arrcube, method, src_satellite, src_sensor, src_LBA,
tgt_satellite, tgt_sensor, tgt_LBA, nodataVal=None, compute_errors=False,
**fallback_argskwargs):
bandwise_errors=True, **fallback_argskwargs):
# type: (Union[np.ndarray, GeoArray], str, str, str, list, str, str, list, int, bool, dict) -> tuple
"""Predict spectral bands of target sensor by applying a machine learning approach.
......@@ -179,6 +180,8 @@ class SpectralHomogenizer(object):
:param nodataVal: no data value
:param compute_errors: whether to compute pixel- / bandwise model errors for estimated pixel values
(default: false)
:param bandwise_errors whether to compute error information for each band separately (True - default)
or to average errors over bands using median (False) (ignored in case of fallback)
:param fallback_argskwargs: arguments and keyword arguments for fallback algorithm ({'args':{}, 'kwargs': {}}
:return: predicted array (rows x columns x bands)
:rtype: Tuple[np.ndarray, Union[np.ndarray, None]]
......@@ -219,11 +222,17 @@ class SpectralHomogenizer(object):
if compute_errors:
errors = PR.compute_prediction_errors(im_homo, cls, nodataVal=nodataVal)
if not bandwise_errors:
errors = np.median(errors, axis=2).astype(errors.dtype)
elif fallback_argskwargs:
# fallback: use linear interpolation and set errors to an array of zeros
im_homo = self.interpolate_cube(arrcube, *fallback_argskwargs['args'], **fallback_argskwargs['kwargs'])
if compute_errors:
errors = np.zeros_like(im_homo, dtype=np.int16)
if not bandwise_errors:
errors = np.zeros_like(im_homo, dtype=np.int16)
else:
errors = np.zeros(im_homo.shape[:2], dtype=np.int16)
else:
raise exc
......
......@@ -19,6 +19,7 @@ dtype_lib_GDAL_Python = {"uint8": 1, "int8": 1, "uint16": 2, "int16": 3, "uint32
proc_chain = ['L1A', 'L1B', 'L1C', 'L2A', 'L2B', 'L2C']
db_jobs_statistics_def = {'pending': 1, 'started': 2, None: 2, 'L1A': 3, 'L1B': 4, 'L1C': 5, 'L2A': 6, 'L2B': 7,
'L2C': 8, 'FAILED': 9} # NOTE: OrderedDicts passed to L1A_map have proc_level=None
bandslist_all_errors = ['ac_errors', 'mask_clouds_confidence', 'spec_homo_errors']
def get_GMS_sensorcode(GMS_identifier):
......
......@@ -45,6 +45,7 @@ class Dataset(object):
self._ac_options = {}
self._ac_errors = None
self._spec_homo_errors = None
self._accuracy_layers = None
# defaults
self.proc_level = 'init'
......@@ -296,32 +297,6 @@ class Dataset(object):
def mask_clouds(self):
self._mask_clouds = None
@property
def mask_clouds_confidence(self):
return self._mask_clouds_confidence
@mask_clouds_confidence.setter
def mask_clouds_confidence(self, *geoArr_initArgs):
if geoArr_initArgs[0] is not None:
cnfArr = GeoArray(*geoArr_initArgs)
assert cnfArr.shape == self.arr.shape[:2], \
"The 'mask_clouds_confidence' GeoArray can only be instanced with an array of the same dimensions " \
"like GMS_obj.arr. Got %s." % str(cnfArr.shape)
if cnfArr._nodata is None:
cnfArr.nodata = DEF_D.get_outFillZeroSaturated(cnfArr.dtype)[0]
cnfArr.gt = self.arr.gt
cnfArr.prj = self.arr.prj
self._mask_clouds_confidence = cnfArr
else:
del self._mask_clouds_confidence
@mask_clouds_confidence.deleter
def mask_clouds_confidence(self):
self._mask_clouds_confidence = None
@property
def dem(self):
"""
......@@ -383,69 +358,6 @@ class Dataset(object):
def pathGen(self):
self._pathGen = None
@property
def ac_errors(self):
"""Returns an instance of GeoArray containing error information calculated by the atmospheric correction.
:return:
"""
return self._ac_errors # FIXME should give a warning if None
@ac_errors.setter
def ac_errors(self, *geoArr_initArgs):
if geoArr_initArgs[0] is not None:
errArr = GeoArray(*geoArr_initArgs)
assert errArr.shape == self.arr.shape, "The 'ac_errors' GeoArray can only be instanced with an array of " \
"the same dimensions like GMS_obj.arr. Got %s." % str(errArr.shape)
if errArr._nodata is None:
errArr.nodata = DEF_D.get_outFillZeroSaturated(errArr.dtype)[0]
errArr.gt = self.arr.gt
errArr.prj = self.arr.prj
errArr.bandnames = self.LBA2bandnames(self.LayerBandsAssignment)
self._ac_errors = errArr
else:
del self.ac_errors
@ac_errors.deleter
def ac_errors(self):
self._ac_errors = None
@property
def spec_homo_errors(self):
"""Returns an instance of GeoArray containing error information calculated during spectral homogenization.
:return:
"""
return self._spec_homo_errors # FIXME should give a warning if None
@spec_homo_errors.setter
def spec_homo_errors(self, *geoArr_initArgs):
if geoArr_initArgs[0] is not None:
errArr = GeoArray(*geoArr_initArgs)
assert errArr.shape == self.arr.shape, "The 'spec_homo_errors' GeoArray can only be instanced with an " \
"array of the same dimensions like GMS_obj.arr. Got %s." \
% str(errArr.shape)
if errArr._nodata is None:
errArr.nodata = DEF_D.get_outFillZeroSaturated(errArr.dtype)[0]
errArr.gt = self.arr.gt
errArr.prj = self.arr.prj
errArr.bandnames = self.LBA2bandnames(self.LayerBandsAssignment)
self._spec_homo_errors = errArr
else:
del self.spec_homo_errors
@spec_homo_errors.deleter
def spec_homo_errors(self):
self._spec_homo_errors = None
@property
def subset(self):
return [self.arr_shape, self.arr_pos]
......@@ -512,15 +424,16 @@ class Dataset(object):
# copy object
sub_GMS_obj = HLP_F.parentObjDict[self.proc_level](*HLP_F.initArgsDict[self.proc_level]) # init
sub_GMS_obj.__dict__ = {k: getattr(self, k) for k in self.__dict__.keys()
if not isinstance(getattr(self, k), (GeoArray, np.ndarray))}.copy()
sub_GMS_obj.__dict__.update(
{k: getattr(self, k) for k in self.__dict__.keys()
if not isinstance(getattr(self, k), (GeoArray, np.ndarray))}.copy())
sub_GMS_obj = copy.deepcopy(sub_GMS_obj)
# clip all array attributes using the given bounds
# list_arraynames = [i for i in self.__dict__ if not callable(getattr(self, i)) and \
# isinstance(getattr(self, i), np.ndarray)]
list_arraynames = [i for i in ['arr', 'masks', 'ac_errors', 'mask_clouds_confidence']
list_arraynames = [i for i in ['arr', 'masks', 'ac_errors', 'mask_clouds_confidence', 'spec_homo_errors']
if hasattr(self, i) and getattr(self, i) is not None] # FIXME hardcoded
assert list_arraynames
assert imBounds or mapBounds, "Either 'imBounds' or 'mapBounds' must be given. Got nothing."
......
......@@ -331,6 +331,164 @@ class GMS_object(Dataset):
def masks(self):
self._masks = None
@property
def mask_clouds_confidence(self):
return self._mask_clouds_confidence
@mask_clouds_confidence.setter
def mask_clouds_confidence(self, *geoArr_initArgs):
if geoArr_initArgs[0] is not None:
cnfArr = GeoArray(*geoArr_initArgs)
assert cnfArr.shape == self.arr.shape[:2], \
"The 'mask_clouds_confidence' GeoArray can only be instanced with an array of the same dimensions " \
"like GMS_obj.arr. Got %s." % str(cnfArr.shape)
if cnfArr._nodata is None:
cnfArr.nodata = DEF_D.get_outFillZeroSaturated(cnfArr.dtype)[0]
cnfArr.gt = self.arr.gt
cnfArr.prj = self.arr.prj
cnfArr.bandnames = ['Confidence of cloud mask']
self._mask_clouds_confidence = cnfArr
else:
del self._mask_clouds_confidence
@mask_clouds_confidence.deleter
def mask_clouds_confidence(self):
self._mask_clouds_confidence = None
@property
def ac_errors(self):
"""Returns an instance of GeoArray containing error information calculated by the atmospheric correction.
:return:
"""
return self._ac_errors # FIXME should give a warning if None
@ac_errors.setter
def ac_errors(self, *geoArr_initArgs):
if geoArr_initArgs[0] is not None:
errArr = GeoArray(*geoArr_initArgs)
if CFG.ac_bandwise_accuracy:
assert errArr.shape == self.arr.shape, \
"The 'ac_errors' GeoArray can only be instanced with an array of the same dimensions like " \
"GMS_obj.arr. Got %s." % str(errArr.shape)
else:
assert errArr.shape[:2] == self.arr.shape[:2], \
"The 'ac_errors' GeoArray can only be instanced with an array of the same X/Y dimensions like " \
"GMS_obj.arr. Got %s." % str(errArr.shape)
if errArr._nodata is None:
errArr.nodata = DEF_D.get_outFillZeroSaturated(errArr.dtype)[0]
errArr.gt = self.arr.gt
errArr.prj = self.arr.prj
errArr.bandnames = \
['AC errors %s' % bN for bN in self.LBA2bandnames(self.LayerBandsAssignment)]\
if CFG.ac_bandwise_accuracy else ['median AC errors']
self._ac_errors = errArr
else:
del self.ac_errors
@ac_errors.deleter
def ac_errors(self):
self._ac_errors = None
@property
def spec_homo_errors(self):
"""Returns an instance of GeoArray containing error information calculated during spectral homogenization.
:return:
"""
return self._spec_homo_errors # FIXME should give a warning if None
@spec_homo_errors.setter
def spec_homo_errors(self, *geoArr_initArgs):
if geoArr_initArgs[0] is not None:
errArr = GeoArray(*geoArr_initArgs)
if CFG.spechomo_bandwise_accuracy:
assert errArr.shape == self.arr.shape, \
"The 'spec_homo_errors' GeoArray can only be instanced with an array of the same dimensions like " \
"GMS_obj.arr. Got %s." % str(errArr.shape)
else:
assert errArr.shape[:2] == self.arr.shape[:2], \
"The 'spec_homo_errors' GeoArray can only be instanced with an array of the same X/Y dimensions " \
"like GMS_obj.arr. Got %s." % str(errArr.shape)
if errArr._nodata is None:
errArr.nodata = DEF_D.get_outFillZeroSaturated(errArr.dtype)[0]
errArr.gt = self.arr.gt
errArr.prj = self.arr.prj
errArr.bandnames = \
['SpecHomo errors %s' % bN for bN in self.LBA2bandnames(self.LayerBandsAssignment)] \
if CFG.spechomo_bandwise_accuracy else ['median SpecHomo errors']
self._spec_homo_errors = errArr
else:
del self.spec_homo_errors
@spec_homo_errors.deleter
def spec_homo_errors(self):
self._spec_homo_errors = None
@property
def accuracy_layers(self):
errs = OrderedDict((band, getattr(self, band)) for band in DEF_D.bandslist_all_errors)
errs_data = OrderedDict((band, err) for band, err in errs.items() if isinstance(err, (np.ndarray, GeoArray)))
bandNames = list(chain.from_iterable([err_data.bandnames for err_data in errs_data.values()]))
if errs_data:
if self._accuracy_layers is not None and list(self._accuracy_layers.bandnames.keys()) == bandNames:
return self._accuracy_layers
else:
self.logger.info('Generating combined accuracy layers array..')
# set accuracy_layers
self._accuracy_layers = GeoArray(np.array(np.dstack(errs_data.values()), dtype='int16'),
geotransform=list(errs_data.values())[0].geotransform,
projection=list(errs_data.values())[0].projection,
bandnames=bandNames,
nodata=DEF_D.get_outFillZeroSaturated('int16')[0])
return self._accuracy_layers
@accuracy_layers.setter
def accuracy_layers(self, geoArr_initArgs):
if geoArr_initArgs[0] is not None:
acc_lay = GeoArray(geoArr_initArgs)
assert acc_lay.shape[:2] == self.arr.shape[:2],\
"The 'accuracy_layers' GeoArray can only be instanced with an array of the same dimensions like " \
"GMS_obj.arr. Got %s." % str(acc_lay.shape)
if acc_lay._nodata is None:
acc_lay.nodata = DEF_D.get_outFillZeroSaturated(acc_lay.dtype)[0]
acc_lay.gt = self.arr.gt
acc_lay.prj = self.arr.prj
if not acc_lay.bandnames:
raise ValueError
self._accuracy_layers = acc_lay
else:
del self._accuracy_layers
@property
def accuracy_layers_meta(self):
if self._accuracy_layers is not None:
return {'map info': geotransform2mapinfo(self._accuracy_layers.gt, self._accuracy_layers.projection),
'coordinate system string': self._accuracy_layers.projection,
'bands': self._accuracy_layers.bands,
'band names': list(self._accuracy_layers.bandnames),
'data ignore value': self._accuracy_layers.nodata}
else:
return None
@property
def cloud_masking_algorithm(self):
if not self._cloud_masking_algorithm:
......@@ -1109,26 +1267,32 @@ class GMS_object(Dataset):
MGRS_info=MGRS_info).get_path_imagedata()
# set dicts
image_type_dict = {'arr': self.image_type, 'masks': 'MAS', 'mask_clouds': 'MAC'}
metaAttr_dict = {'arr': 'meta_odict', 'masks': 'masks_meta', 'mask_clouds': 'masks_meta'}
out_dtype_dict = {'arr': 'int16', 'masks': 'uint8', 'mask_clouds': 'uint8'}
print_dict = {'RSD_L1A': 'L1A satellite data', 'MAS_L1A': 'L1A masks', 'MAC_L1A': 'L1A cloud mask',
'RSD_L1B': 'L1B satellite data', 'MAS_L1B': 'L1B masks', 'MAC_L1B': 'L1B cloud mask',
'RSD_L1C': 'L1C atm. corrected reflectance data', 'MAS_L1C': 'L1C masks',
'MAC_L1C': 'L1C cloud mask',
'RSD_L2A': 'L2A geometrically homogenized data', 'MAS_L2A': 'L2A masks',
'MAC_L2A': 'L2A cloud mask',
'RSD_L2B': 'L2B spectrally homogenized data', 'MAS_L2B': 'L2B masks', 'MAC_L2B': 'L2B cloud mask',
'RSD_L2C': 'L2C MGRS tiled data', 'MAS_L2C': 'L2C masks', 'MAC_L2C': 'L2C cloud mask'}
image_type_dict = dict(arr=self.image_type, masks='MAS', mask_clouds='MAC', accuracy_layers='AL')
metaAttr_dict = dict(arr='meta_odict', masks='masks_meta', mask_clouds='masks_meta',
accuracy_layers='accuracy_layers_meta')
out_dtype_dict = dict(arr='int16', masks='uint8', mask_clouds='uint8', accuracy_layers='int16')
print_dict = dict(
RSD_L1A='satellite data', MAS_L1A='masks', MAC_L1A='cloud mask',
RSD_L1B='satellite data', MAS_L1B='masks', MAC_L1B='cloud mask',
RSD_L1C='atm. corrected reflectance data', MAS_L1C='masks', MAC_L1C='cloud mask', AL_L1C='accuracy layers',
RSD_L2A='geometrically homogenized data', MAS_L2A='masks', MAC_L2A='cloud mask', AL_L2A='accuracy layers',
RSD_L2B='spectrally homogenized data', MAS_L2B='masks', MAC_L2B='cloud mask', AL_L2B='accuracy layers',
RSD_L2C='MGRS tiled data', MAS_L2C='masks', MAC_L2C='cloud mask', AL_L2C='accuracy layers',)
# ensure self.masks exists
if not hasattr(self, 'masks') or self.masks is None:
self.build_combined_masks_array() # creates InObj.masks and InObj.masks_meta
# loop through all attributes to write and execute writer
attributes2write = ['arr', 'masks'] + (['mask_clouds'] if write_masks_as_ENVI_classification else [])
attributes2write = attributes2write if self.proc_level not in ['L1A', 'L1B'] else ['arr', 'masks']
# generate list of attributes to write
attributes2write = ['arr', 'masks']
if self.proc_level not in ['L1A', 'L1B'] and write_masks_as_ENVI_classification:
attributes2write.append('mask_clouds')
if CFG.ac_estimate_accuracy or CFG.spechomo_estimate_accuracy:
attributes2write.append('accuracy_layers')
# loop through all attributes to write and execute writer
for arrayname in attributes2write:
descriptor = '%s_%s' % (image_type_dict[arrayname], self.proc_level)
......@@ -1149,7 +1313,7 @@ class GMS_object(Dataset):
meta_attr = metaAttr_dict[arrayname]
if not is_tempfile:
self.log_for_fullArr_or_firstTile('Writing %s.' % print_dict[descriptor])
self.log_for_fullArr_or_firstTile('Writing %s %s.' % (self.proc_level, print_dict[descriptor]))
#########################
# GeoArray in disk mode #
......
......@@ -362,6 +362,8 @@ class JobConfig(object):
gp('ac_max_ram_gb', json_processors['L1C']['ac_max_ram_gb'])
self.ac_estimate_accuracy = \
gp('ac_estimate_accuracy', json_processors['L1C']['ac_estimate_accuracy'])
self.ac_bandwise_accuracy = \
gp('ac_bandwise_accuracy', json_processors['L1C']['ac_bandwise_accuracy'])
# L2A
self.exec_L2AP = gp('exec_L2AP', [
......@@ -381,6 +383,8 @@ class JobConfig(object):
self.spechomo_method = gp('spechomo_method', json_processors['L2B']['spechomo_method'])
self.spechomo_estimate_accuracy = \
gp('spechomo_estimate_accuracy', json_processors['L2B']['spechomo_estimate_accuracy'])
self.spechomo_bandwise_accuracy = \
gp('spechomo_bandwise_accuracy', json_processors['L2B']['spechomo_bandwise_accuracy'])
# L2C
self.exec_L2CP = gp('exec_L2CP', [
......
......@@ -100,7 +100,9 @@
["Clear", "Snow", "Water", "Shadow", "Cirrus", "Cloud"] => whole image is atmosperically corrected*/
"ac_scale_factor_errors": 255,
"ac_max_ram_gb": 20, /*maximum amount of RAM to be allocated for atmospheric correction [gigabytes]*/
"ac_estimate_accuracy": false /*whether to produce an 'AC errors' and a 'mask confidence' array*/
"ac_estimate_accuracy": false, /*whether to produce an 'AC errors' and a 'mask confidence' array*/
"ac_bandwise_accuracy": false /*if true, the 'AC_errors' array is produced bandwise, otherwise bands are
averaged using median*/
},
"L2A": { /*Level 2A processing: geometric homogenization*/
......@@ -121,8 +123,10 @@
"delete_output": false,
"spechomo_method": "LR", /*Method used for spectral homogenization.
/*LI: Linear interpolation; LR: Linear regression; RR: Ridge regression*/
"spechomo_estimate_accuracy": false /*whether to produce pixel- and bandwise information about estimation
acurracy of spectral homogenization*/
"spechomo_estimate_accuracy": false, /*whether to produce pixel- and bandwise information about estimation
acurracy of spectral homogenization*/
"spechomo_bandwise_accuracy": true /*if true, the spectral homogenization accuracy array is produced
bandwise, otherwise accuracy information is averaged using median*/
},
"L2C": {
......
......@@ -97,6 +97,7 @@ gms_schema_input = dict(
ac_scale_factor_errors=dict(type='integer', required=False),
ac_max_ram_gb=dict(type='integer', required=False),
ac_estimate_accuracy=dict(type='boolean', required=False),
ac_bandwise_accuracy=dict(type='boolean', required=False),
)),
L2A=dict(type='dict', required=False, schema=dict(
run_processor=dict(type='boolean', required=False),
......@@ -115,6 +116,7 @@ gms_schema_input = dict(
delete_output=dict(type='boolean', required=False),
spechomo_method=dict(type='string', required=False, allowed=['LI', 'LR', 'RR']),
spechomo_estimate_accuracy=dict(type='boolean', required=False),
spechomo_bandwise_accuracy=dict(type='boolean', required=False),
)),
L2C=dict(type='dict', required=False, schema=dict(
run_processor=dict(type='boolean', required=False),
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment