Commit 726e5b1a authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Second version of accuracy layers + writer. Accuracy layers are now only written during L2C_P.

Added L2C_P.AccuracyCube().
Added logging to DEM_Creator.
parent 37a7f8c9
......@@ -57,6 +57,7 @@ class L2B_object(L2A_object):
#################################################################
# collect some information specifying the needed homogenization #
#################################################################
method = CFG.spechomo_method
src_dsID = sat_sen_to_datasetid(self.satellite, self.sensor)
src_cwls = self.meta_odict['wavelength']
......@@ -71,6 +72,7 @@ class L2B_object(L2A_object):
####################################################
# special cases where homogenization is not needed #
####################################################
if self.dataset_ID == CFG.datasetid_spectral_ref:
self.logger.info("Spectral homogenization has been skipped because the dataset id equals the dataset id of "
"the spectral refernce sensor.")
......@@ -82,15 +84,22 @@ class L2B_object(L2A_object):
"are already equal to the target sensor's.")
return
####################################
# perform spectral homogenization! #
####################################
#################################################
# perform spectral homogenization of image data #
#################################################
SpH = SpectralHomogenizer(classifier_rootDir=CFG.path_spechomo_classif, logger=self.logger)
if method == 'LI' or CFG.datasetid_spectral_ref is None:
# linear interpolation (if intended by user or in case of custom spectral characteristics of target sensor)
# -> no classifier for that case available -> linear interpolation
im, errs = SpH.interpolate_cube(self.arr, src_cwls, tgt_cwls, kind='linear')
im = SpH.interpolate_cube(self.arr, src_cwls, tgt_cwls, kind='linear')
if CFG.spechomo_estimate_accuracy:
self.logger.warning("Unable to compute any error information in case spectral homogenization algorithm "
"is set to 'LI' (Linear Interpolation)")
errs = None
else:
# a known sensor has been specified as spectral reference => apply a machine learner
......@@ -113,6 +122,7 @@ class L2B_object(L2A_object):
###################
# update metadata #
###################
self.LayerBandsAssignment = tgt_LBA
self.meta_odict['wavelength'] = list(tgt_cwls)
self.meta_odict['bands'] = len(tgt_cwls)
......@@ -122,6 +132,17 @@ class L2B_object(L2A_object):
self.arr = im # type: GeoArray
self.spec_homo_errors = errs # type: Union[np.ndarray, None] # int16, None if ~CFG.spechomo_estimate_accuracy
#########################################################################################
# perform spectral homogenization of bandwise error information from earlier processors #
#########################################################################################
if self.ac_errors and self.ac_errors.ndim == 3:
self.logger.info("Performing linear interpolation for 'AC errors' array to match target sensor bands "
"number..")
outarr = \
interp1d(np.array(src_cwls), self.ac_errors, axis=2, kind='linear', fill_value='extrapolate')(tgt_cwls)
self.ac_errors = outarr.astype(np.int16)
class SpectralHomogenizer(object):
"""Class for applying spectral homogenization by applying an interpolation or machine learning approach."""
......@@ -229,6 +250,8 @@ class SpectralHomogenizer(object):
# 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:
self.logger.warning("Spectral homogenization algorithm had to be performed by linear interpolation "
"(fallback). Unable to compute any accuracy information from that.")
if not bandwise_errors:
errors = np.zeros_like(im_homo, dtype=np.int16)
else:
......
# -*- coding: utf-8 -*-
"""Level 2C Processor: Quality layers"""
import numpy as np
from collections import OrderedDict
from geoarray import GeoArray
from ..misc.definition_dicts import bandslist_all_errors, get_outFillZeroSaturated
from .L2B_P import L2B_object
from ..options.config import GMS_config as CFG
__author__ = 'Daniel Scheffler'
......@@ -20,8 +27,71 @@ class L2C_object(L2B_object):
self.proc_level = 'L2C'
self.proc_status = 'initialized'
def calc_geometric_accurracy(self):
pass
def calc_spectral_accurracy(self):
pass
class AccuracyCube(GeoArray):
def __init__(self, GMS_obj):
self._GMS_obj = GMS_obj
# privates
self._layers = None
if self.layers:
super(AccuracyCube, self).__init__(self.generate_array(),
geotransform=list(self.layers.values())[0].gt,
projection=list(self.layers.values())[0].prj,
bandnames=self.get_bandnames(),
nodata=get_outFillZeroSaturated('int16')[0])
else:
# return an empty GeoArray
super(AccuracyCube, self).__init__(np.array([]))
@property
def layers(self):
# type: () -> OrderedDict
if not self._layers:
errs = OrderedDict((band, getattr(self._GMS_obj, band)) for band in bandslist_all_errors)
self._layers = \
OrderedDict((band, err) for band, err in errs.items() if isinstance(err, (np.ndarray, GeoArray)))
return self._layers
def get_bandnames(self):
bandnames = []
for errArrName, errArr in self.layers.items():
if errArrName == 'ac_errors':
if CFG.ac_bandwise_accuracy:
bandnames.extend(['AC errors %s' % bN for bN in errArr.bandnames])
else:
bandnames.append('median of AC errors')
elif errArrName == 'mask_clouds_confidence':
bandnames.append('confidence of cloud mask')
elif errArrName == 'spec_homo_errors':
if CFG.spechomo_bandwise_accuracy:
bandnames.extend(['error of spectral homogenization %s' % bN for bN in errArr.bandnames])
else:
bandnames.append('median error of spectral homogenization')
else:
raise RuntimeError('Error setting bandnames for %s.' % errArrName)
return bandnames
def generate_array(self):
err_layers = self.layers.copy() # copy OrdDict, otherwise attributes of GMS_object are overwritten
# handle CFG.ac_bandwise_accuracy (average ac_errors if needed)
if 'ac_errors' in err_layers and not CFG.ac_bandwise_accuracy and err_layers['ac_errors'].bands > 1:
err_layers['ac_errors'].arr = \
np.median(err_layers['ac_errors'], axis=2).astype(err_layers['ac_errors'].dtype)
# handle CFG.spechomo_bandwise_accuracy (average spec_homo_errors if needed)
if 'spec_homo_errors' in err_layers and not CFG.spechomo_bandwise_accuracy and \
err_layers['spec_homo_errors'].bands > 1:
err_layers['spec_homo_errors'].arr = \
np.median(err_layers['spec_homo_errors'], axis=2).astype(err_layers['spec_homo_errors'].dtype)
return np.array(np.dstack(err_layers.values()), dtype='int16')
......@@ -14,6 +14,7 @@ from logging import Logger
from matplotlib import pyplot as plt
from typing import Union, Dict, List, Tuple # noqa F401 # flake8 issue
from datetime import datetime
import logging
import dill
import numpy as np
......@@ -446,7 +447,7 @@ def open_specific_file_within_archive(path_archive, matching_expression, read_mo
class DEM_Creator(object):
def __init__(self, dem_sensor='SRTM', db_conn=''):
def __init__(self, dem_sensor='SRTM', db_conn='', logger=None):
"""Creator class for digital elevation models based on ASTER or SRTM.
:param dem_sensor: 'SRTM' or 'ASTER'
......@@ -457,6 +458,7 @@ class DEM_Creator(object):
self.dem_sensor = dem_sensor
self.db_conn = db_conn if db_conn else CFG.conn_database
self.logger = logger or logging.getLogger('DEM_Creator')
self.project_dir = os.path.abspath(os.path.curdir)
self.rootpath_DEMs = ''
......@@ -521,10 +523,13 @@ class DEM_Creator(object):
sceneIDs_srtm = self.get_overlapping_DEM_tiles(tgt_corner_coord_lonlat, timeout_sec)
except ConnectionRefusedError:
# fallback to plain pgSQL
self.logger.warning('SpatialIndexMediator refused connection. Falling back to plain postgreSQL query.')
sceneIDs_srtm = self.get_overlapping_DEM_tiles(tgt_corner_coord_lonlat, use_index_mediator=False)
if not sceneIDs_srtm:
# fallback to plain pgSQL
self.logger.warning('SpatialIndexMediator did not return matching DEM tiles. '
'Trying plain postgreSQL query..')
sceneIDs_srtm = self.get_overlapping_DEM_tiles(tgt_corner_coord_lonlat, use_index_mediator=False)
except TimeoutError:
......@@ -544,13 +549,12 @@ class DEM_Creator(object):
return gdalImPaths
@staticmethod
def _run_cmd(cmd):
def _run_cmd(self, cmd):
output, exitcode, err = HLP_F.subcall_with_output(cmd)
if exitcode:
print('\nAn error occurred during DEM creation.')
print("Print traceback in case you care:")
print(err.decode('latin-1'))
self.logger.error('\nAn error occurred during DEM creation.')
self.logger.warning("Print traceback in case you care:")
self.logger.warning(err.decode('latin-1'))
if output:
return output.decode('UTF-8')
......
......@@ -307,11 +307,11 @@ class Dataset(object):
self.logger.info('Generating DEM...')
DC_args = (self.arr.box.boxMapXY, self.arr.prj, self.arr.xgsd, self.arr.ygsd)
try:
DC = INP_R.DEM_Creator(dem_sensor='SRTM')
DC = INP_R.DEM_Creator(dem_sensor='SRTM', logger=self.logger)
self._dem = DC.from_extent(*DC_args)
except RuntimeError:
self.logger.warning('SRTM DEM generation failed. Falling back to ASTER...')
DC = INP_R.DEM_Creator(dem_sensor='ASTER')
DC = INP_R.DEM_Creator(dem_sensor='ASTER', logger=self.logger)
self._dem = DC.from_extent(*DC_args)
self._dem.nodata = DEF_D.get_outFillZeroSaturated(self._dem.dtype)[0]
......
......@@ -348,7 +348,7 @@ class GMS_object(Dataset):
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']
cnfArr.bandnames = ['confidence']
self._mask_clouds_confidence = cnfArr
else:
......@@ -385,9 +385,7 @@ class GMS_object(Dataset):
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']
errArr.bandnames = self.LBA2bandnames(self.LayerBandsAssignment) if errArr.ndim == 3 else ['median']
self._ac_errors = errArr
else:
......@@ -424,9 +422,7 @@ class GMS_object(Dataset):
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']
errArr.bandnames = self.LBA2bandnames(self.LayerBandsAssignment) if errArr.ndim == 3 else ['median']
self._spec_homo_errors = errArr
else:
......@@ -438,32 +434,19 @@ class GMS_object(Dataset):
@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 not self._accuracy_layers:
if not self.proc_level.startswith('L2'):
self.logger.warning('Attempt to get %s accuracy layers failed - they are a Level 2 feature only.'
% self.proc_level)
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..')
if 'ac_errors' in errs_data and self.proc_level == 'L1C' and not CFG.ac_bandwise_accuracy:
errs_data['ac_errors'] = GeoArray(np.median(errs_data['ac_errors'], axis=2),
geotransform=errs_data['ac_errors'].geotransform,
projection=errs_data['ac_errors'].projection,
bandnames=['median AC errors']
)
# 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])
self.logger.info('Generating combined accuracy layers array..')
try:
from ..algorithms.L2C_P import AccuracyCube
self._accuracy_layers = AccuracyCube(self)
except Exception as e:
raise RuntimeError('Failed to generate AccuracyCube!', e)
return self._accuracy_layers
return self._accuracy_layers
@accuracy_layers.setter
def accuracy_layers(self, geoArr_initArgs):
......@@ -485,6 +468,10 @@ class GMS_object(Dataset):
else:
del self._accuracy_layers
@accuracy_layers.deleter
def accuracy_layers(self):
self._accuracy_layers = None
@property
def accuracy_layers_meta(self):
if self._accuracy_layers is not None:
......@@ -1283,9 +1270,9 @@ class GMS_object(Dataset):
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_L1C='atm. corrected reflectance data', MAS_L1C='masks', MAC_L1C='cloud mask',
RSD_L2A='geometrically homogenized data', MAS_L2A='masks', MAC_L2A='cloud mask',
RSD_L2B='spectrally homogenized data', MAS_L2B='masks', MAC_L2B='cloud mask',
RSD_L2C='MGRS tiled data', MAS_L2C='masks', MAC_L2C='cloud mask', AL_L2C='accuracy layers',)
# ensure self.masks exists
......@@ -1298,10 +1285,13 @@ class GMS_object(Dataset):
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')
if self.proc_level == 'L2C': # and CFG.ac_estimate_accuracy or CFG.spechomo_estimate_accuracy:
attributes2write.append('accuracy_layers')
###########################################################
# loop through all attributes to write and execute writer #
###########################################################
# loop through all attributes to write and execute writer
for arrayname in attributes2write:
descriptor = '%s_%s' % (image_type_dict[arrayname], self.proc_level)
......@@ -1485,6 +1475,7 @@ class GMS_object(Dataset):
outpath_arr = os.path.splitext(outpath_hdr)[0] + '.%s' % self.outInterleave
if not CFG.inmem_serialization:
setattr(self, arrayname, outpath_arr) # replace array by output path
if arrayname == 'masks':
setattr(self, 'mask_nodata', outpath_arr)
......@@ -1640,6 +1631,7 @@ class GMS_object(Dataset):
del self.masks
del self.mask_clouds_confidence
del self.ac_errors
del self.accuracy_layers
class failed_GMS_object(GMS_object):
......
......@@ -124,9 +124,11 @@
"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*/
acurracy of spectral homogenization
NOTE: only possible if 'spechomo_method' is not linear interpol.*/
"spechomo_bandwise_accuracy": true /*if true, the spectral homogenization accuracy array is produced
bandwise, otherwise accuracy information is averaged using median*/
},
"L2C": {
......
......@@ -517,7 +517,7 @@ class process_controller(object):
# clear any temporary files
tempdir = os.path.join(self.config.path_tempdir)
self.logger.warning('Deleting temporary directory %s.' % tempdir)
self.logger.info('Deleting temporary directory %s.' % tempdir)
if os.path.exists(tempdir):
shutil.rmtree(tempdir)
......
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