Commit 357719f5 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Merge branch 'master' into feature/add_locks

parents e475587d 31828c1a
...@@ -40,6 +40,8 @@ class L1A_object(GMS_object): ...@@ -40,6 +40,8 @@ class L1A_object(GMS_object):
""" """
# TODO docstring # TODO docstring
# NOTE: kwargs is in here to allow instanciating with additional arg 'proc_level'
# load default attribute values and methods # load default attribute values and methods
super(L1A_object, self).__init__() super(L1A_object, self).__init__()
...@@ -274,7 +276,6 @@ class L1A_object(GMS_object): ...@@ -274,7 +276,6 @@ class L1A_object(GMS_object):
Works for: RapidEye (metadata.xml),SPOT(metadata.dim),LANDSAT(mtl.txt),ASTER(downloaded coremetadata), Works for: RapidEye (metadata.xml),SPOT(metadata.dim),LANDSAT(mtl.txt),ASTER(downloaded coremetadata),
ALOS(summary.txt & Leader file) ALOS(summary.txt & Leader file)
:param v:
:return: :return:
""" """
...@@ -522,7 +523,6 @@ class L1A_object(GMS_object): ...@@ -522,7 +523,6 @@ class L1A_object(GMS_object):
if rasObj.get_projection_type() == 'UTM': if rasObj.get_projection_type() == 'UTM':
self.MetaObj.CornerTieP_UTM = rasObj.get_corner_coordinates('UTM') self.MetaObj.CornerTieP_UTM = rasObj.get_corner_coordinates('UTM')
self.meta_odict = self.MetaObj.to_odict() # important in order to keep geotransform/projection
if CFG.inmem_serialization: if CFG.inmem_serialization:
self.delete_tempFiles() # these files are needed later in Python execution mode self.delete_tempFiles() # these files are needed later in Python execution mode
self.MetaObj.Dataname = previous_dataname # /vsi.. pointing directly to raw data archive (which exists) self.MetaObj.Dataname = previous_dataname # /vsi.. pointing directly to raw data archive (which exists)
......
...@@ -16,7 +16,8 @@ import numpy as np ...@@ -16,7 +16,8 @@ import numpy as np
from geopandas import GeoDataFrame from geopandas import GeoDataFrame
from shapely.geometry import box from shapely.geometry import box
import pytz import pytz
from typing import Union # noqa F401 # flake8 issue import traceback
from typing import Union, TYPE_CHECKING # noqa F401 # flake8 issue
from arosics import COREG, DESHIFTER from arosics import COREG, DESHIFTER
from geoarray import GeoArray from geoarray import GeoArray
...@@ -36,6 +37,10 @@ from ..misc.logging import GMS_logger ...@@ -36,6 +37,10 @@ from ..misc.logging import GMS_logger
from ..misc.spatial_index_mediator import SpatialIndexMediator from ..misc.spatial_index_mediator import SpatialIndexMediator
from ..misc.definition_dicts import get_GMS_sensorcode, get_outFillZeroSaturated from ..misc.definition_dicts import get_GMS_sensorcode, get_outFillZeroSaturated
if TYPE_CHECKING:
from shapely.geometry import Polygon # noqa F401 # flake8 issue
from logging import Logger # noqa F401 # flake8 issue
__author__ = 'Daniel Scheffler' __author__ = 'Daniel Scheffler'
...@@ -44,6 +49,7 @@ class Scene_finder(object): ...@@ -44,6 +49,7 @@ class Scene_finder(object):
def __init__(self, src_boundsLonLat, src_AcqDate, src_prj, src_footprint_poly, sceneID_excluded=None, def __init__(self, src_boundsLonLat, src_AcqDate, src_prj, src_footprint_poly, sceneID_excluded=None,
min_overlap=20, min_cloudcov=0, max_cloudcov=20, plusminus_days=30, plusminus_years=10, logger=None): min_overlap=20, min_cloudcov=0, max_cloudcov=20, plusminus_days=30, plusminus_years=10, logger=None):
# type: (list, datetime, str, Polygon, int, int, int, int, int, int, Logger) -> None
"""Initialize Scene_finder. """Initialize Scene_finder.
:param src_boundsLonLat: :param src_boundsLonLat:
...@@ -271,7 +277,7 @@ class L1B_object(L1A_object): ...@@ -271,7 +277,7 @@ class L1B_object(L1A_object):
def get_spatial_reference_scene(self): def get_spatial_reference_scene(self):
boundsLonLat = corner_coord_to_minmax(self.trueDataCornerLonLat) boundsLonLat = corner_coord_to_minmax(self.trueDataCornerLonLat)
footprint_poly = HLP_F.CornerLonLat_to_shapelyPoly(self.trueDataCornerLonLat) footprint_poly = HLP_F.CornerLonLat_to_shapelyPoly(self.trueDataCornerLonLat)
RSF = Scene_finder(boundsLonLat, self.acq_datetime, self.meta_odict['coordinate system string'], RSF = Scene_finder(boundsLonLat, self.acq_datetime, self.MetaObj.projection,
footprint_poly, self.scene_ID, footprint_poly, self.scene_ID,
min_overlap=CFG.spatial_ref_min_overlap, min_overlap=CFG.spatial_ref_min_overlap,
min_cloudcov=CFG.spatial_ref_min_cloudcov, min_cloudcov=CFG.spatial_ref_min_cloudcov,
...@@ -480,13 +486,14 @@ class L1B_object(L1A_object): ...@@ -480,13 +486,14 @@ class L1B_object(L1A_object):
ref_obj = GMS_object.from_disk((path_gmsFile, ['cube', None])) ref_obj = GMS_object.from_disk((path_gmsFile, ['cube', None]))
# get spectral characteristics # get spectral characteristics
ref_cwl, shift_cwl = [[float(i) for i in GMS_obj.meta_odict['wavelength']] for GMS_obj in [ref_obj, self]] ref_cwl = [float(ref_obj.MetaObj.CWL[bN]) for bN in ref_obj.MetaObj.LayerBandsAssignment]
ref_fwhm, shift_fwhm = [[float(i) for i in GMS_obj.meta_odict['bandwidths']] for GMS_obj in [ref_obj, self]] shift_cwl = [float(self.MetaObj.CWL[bN]) for bN in self.MetaObj.LayerBandsAssignment]
shift_fwhm = [float(self.MetaObj.FWHM[bN]) for bN in self.MetaObj.LayerBandsAssignment]
# exclude cirrus/oxygen band of Landsat-8/Sentinel-2 # exclude cirrus/oxygen band of Landsat-8/Sentinel-2
shift_bbl, ref_bbl = [False] * len(shift_cwl), [False] * len(ref_cwl) # bad band lists shift_bbl, ref_bbl = [False] * len(shift_cwl), [False] * len(ref_cwl) # bad band lists
for GMS_obj, s_r, bbl in zip([self, ref_obj], ['shift', 'ref'], [shift_bbl, ref_bbl]): for GMS_obj, s_r, bbl in zip([self, ref_obj], ['shift', 'ref'], [shift_bbl, ref_bbl]):
GMS_obj.GMS_identifier['logger'] = None # set a dummy value in order to avoid Exception GMS_obj.GMS_identifier.logger = None # set a dummy value in order to avoid Exception
sensorcode = get_GMS_sensorcode(GMS_obj.GMS_identifier) sensorcode = get_GMS_sensorcode(GMS_obj.GMS_identifier)
if sensorcode in ['LDCM', 'S2A', 'S2B'] and '9' in GMS_obj.LayerBandsAssignment: if sensorcode in ['LDCM', 'S2A', 'S2B'] and '9' in GMS_obj.LayerBandsAssignment:
bbl[GMS_obj.LayerBandsAssignment.index('9')] = True bbl[GMS_obj.LayerBandsAssignment.index('9')] = True
...@@ -541,6 +548,9 @@ class L1B_object(L1A_object): ...@@ -541,6 +548,9 @@ class L1B_object(L1A_object):
self.logger.warning('Coregistration skipped according to user configuration.') self.logger.warning('Coregistration skipped according to user configuration.')
elif self.coreg_needed and self.spatRef_available: elif self.coreg_needed and self.spatRef_available:
self.coreg_info.update({'reference scene ID': self.spatRef_scene.scene_ID})
self.coreg_info.update({'reference entity ID': self.spatRef_scene.entity_ID})
geoArr_ref = GeoArray(self.spatRef_scene.filePath) geoArr_ref = GeoArray(self.spatRef_scene.filePath)
geoArr_shift = GeoArray(self.arr) geoArr_shift = GeoArray(self.arr)
r_b4match, s_b4match = self.get_opt_bands4matching(target_cwlPos_nm=CFG.coreg_band_wavelength_for_matching) r_b4match, s_b4match = self.get_opt_bands4matching(target_cwlPos_nm=CFG.coreg_band_wavelength_for_matching)
...@@ -552,7 +562,7 @@ class L1B_object(L1A_object): ...@@ -552,7 +562,7 @@ class L1B_object(L1A_object):
max_shift=CFG.coreg_max_shift_allowed, max_shift=CFG.coreg_max_shift_allowed,
ws=CFG.coreg_window_size, ws=CFG.coreg_window_size,
data_corners_ref=[[x, y] for x, y in self.spatRef_scene.polyUTM.convex_hull.exterior.coords], data_corners_ref=[[x, y] for x, y in self.spatRef_scene.polyUTM.convex_hull.exterior.coords],
data_corners_tgt=[transform_any_prj(EPSG2WKT(4326), self.meta_odict['coordinate system string'], x, y) data_corners_tgt=[transform_any_prj(EPSG2WKT(4326), self.MetaObj.projection, x, y)
for x, y in HLP_F.reorder_CornerLonLat(self.trueDataCornerLonLat)], for x, y in HLP_F.reorder_CornerLonLat(self.trueDataCornerLonLat)],
nodata=(get_outFillZeroSaturated(geoArr_ref.dtype)[0], nodata=(get_outFillZeroSaturated(geoArr_ref.dtype)[0],
get_outFillZeroSaturated(geoArr_shift.dtype)[0]), get_outFillZeroSaturated(geoArr_shift.dtype)[0]),
...@@ -561,26 +571,36 @@ class L1B_object(L1A_object): ...@@ -561,26 +571,36 @@ class L1B_object(L1A_object):
q=True q=True
) )
COREG_obj = COREG(geoArr_ref, geoArr_shift, **coreg_kwargs) # initialize COREG object
COREG_obj.calculate_spatial_shifts() try:
COREG_obj = COREG(geoArr_ref, geoArr_shift, **coreg_kwargs)
self.coreg_info.update( except Exception as e:
COREG_obj.coreg_info) # no clipping to trueCornerLonLat until here -> only shift correction COREG_obj = None
self.coreg_info.update({'reference scene ID': self.spatRef_scene.scene_ID}) self.logger.error('\nAn error occurred during coregistration. BE AWARE THAT THE SCENE %s '
self.coreg_info.update({'reference entity ID': self.spatRef_scene.entity_ID}) '(ENTITY ID %s) HAS NOT BEEN COREGISTERED! Error message was: \n%s\n'
self.coreg_info.update({'shift_reliability': COREG_obj.shift_reliability}) % (self.scene_ID, self.entity_ID, repr(e)))
self.logger.error(traceback.format_exc())
if COREG_obj.success: # TODO include that in the job summary
self.coreg_info['success'] = True
self.logger.info("Calculated map shifts (X,Y): %s / %s" # calculate_spatial_shifts
% (COREG_obj.x_shift_map, if COREG_obj:
COREG_obj.y_shift_map)) # FIXME direkt in calculate_spatial_shifts loggen COREG_obj.calculate_spatial_shifts()
self.logger.info("Reliability of calculated shift: %.1f percent" % COREG_obj.shift_reliability)
self.coreg_info.update(
COREG_obj.coreg_info) # no clipping to trueCornerLonLat until here -> only shift correction
self.coreg_info.update({'shift_reliability': COREG_obj.shift_reliability})
if COREG_obj.success:
self.coreg_info['success'] = True
self.logger.info("Calculated map shifts (X,Y): %s / %s"
% (COREG_obj.x_shift_map,
COREG_obj.y_shift_map)) # FIXME direkt in calculate_spatial_shifts loggen
self.logger.info("Reliability of calculated shift: %.1f percent" % COREG_obj.shift_reliability)
else: else:
# TODO add database entry with error hint # TODO add database entry with error hint
[self.logger.error('ERROR during coregistration of scene %s (entity ID %s):\n%s' [self.logger.error('ERROR during coregistration of scene %s (entity ID %s):\n%s'
% (self.scene_ID, self.entity_ID, err)) for err in COREG_obj.tracked_errors] % (self.scene_ID, self.entity_ID, err)) for err in COREG_obj.tracked_errors]
else: else:
if self.coreg_needed: if self.coreg_needed:
...@@ -671,10 +691,10 @@ class L1B_object(L1A_object): ...@@ -671,10 +691,10 @@ class L1B_object(L1A_object):
# update geoinformations and array shape related attributes # update geoinformations and array shape related attributes
self.logger.info("Updating geoinformations of '%s' attribute..." % attrname) self.logger.info("Updating geoinformations of '%s' attribute..." % attrname)
if attrname == 'arr': if attrname == 'arr':
self.meta_odict['map info'] = DS.updated_map_info self.MetaObj.map_info = DS.updated_map_info
self.meta_odict['coordinate system string'] = EPSG2WKT(WKT2EPSG(DS.updated_projection)) self.MetaObj.projection = EPSG2WKT(WKT2EPSG(DS.updated_projection))
self.shape_fullArr = DS.arr_shifted.shape self.shape_fullArr = DS.arr_shifted.shape
self.meta_odict['lines'], self.meta_odict['samples'] = DS.arr_shifted.shape[:2] self.MetaObj.rows, self.MetaObj.cols = DS.arr_shifted.shape[:2]
else: else:
self.masks_meta['map info'] = DS.updated_map_info self.masks_meta['map info'] = DS.updated_map_info
self.masks_meta['coordinate system string'] = EPSG2WKT(WKT2EPSG(DS.updated_projection)) self.masks_meta['coordinate system string'] = EPSG2WKT(WKT2EPSG(DS.updated_projection))
...@@ -686,7 +706,7 @@ class L1B_object(L1A_object): ...@@ -686,7 +706,7 @@ class L1B_object(L1A_object):
geoArr.arr, geoArr.gt, geoArr.prj = \ geoArr.arr, geoArr.gt, geoArr.prj = \
DS.GeoArray_shifted.arr, DS.GeoArray_shifted.gt, DS.GeoArray_shifted.prj DS.GeoArray_shifted.arr, DS.GeoArray_shifted.gt, DS.GeoArray_shifted.prj
# setattr(self,attrname, DS.GeoArray_shifted) # NOTE: don't set array earlier because setter will also # setattr(self,attrname, DS.GeoArray_shifted) # NOTE: don't set array earlier because setter will also
# # update arr.gt/.prj/.nodata from meta_odict # # update arr.gt/.prj/.nodata from MetaObj
self.resamp_needed = False self.resamp_needed = False
self.coreg_needed = False self.coreg_needed = False
......
...@@ -64,8 +64,8 @@ class L1C_object(L1B_object): ...@@ -64,8 +64,8 @@ class L1C_object(L1B_object):
self.logger.info('Calculating LonLat array...') self.logger.info('Calculating LonLat array...')
self._lonlat_arr = \ self._lonlat_arr = \
GEOP.get_lonlat_coord_array(self.shape_fullArr, self.arr_pos, GEOP.get_lonlat_coord_array(self.shape_fullArr, self.arr_pos,
mapinfo2geotransform(self.meta_odict['map info']), mapinfo2geotransform(self.MetaObj.map_info),
self.meta_odict['coordinate system string'], self.MetaObj.projection,
meshwidth=10, # for faster processing meshwidth=10, # for faster processing
nodata_mask=None, # dont overwrite areas outside the image with nodata nodata_mask=None, # dont overwrite areas outside the image with nodata
outFill=get_outFillZeroSaturated(np.float32)[0])[0] outFill=get_outFillZeroSaturated(np.float32)[0])[0]
...@@ -87,17 +87,18 @@ class L1C_object(L1B_object): ...@@ -87,17 +87,18 @@ class L1C_object(L1B_object):
""" """
if self._VZA_arr is None: if self._VZA_arr is None:
self.logger.info('Calculating viewing zenith array...') self.logger.info('Calculating viewing zenith array...')
if 'ViewingAngle_arrProv' in self.meta_odict and self.meta_odict['ViewingAngle_arrProv']: if self.MetaObj.ViewingAngle_arrProv:
# Sentinel-2 # Sentinel-2
self._VZA_arr = GEOP.adjust_acquisArrProv_to_shapeFullArr(self.meta_odict['ViewingAngle_arrProv'], self._VZA_arr = GEOP.adjust_acquisArrProv_to_shapeFullArr(
self.shape_fullArr, {k: v.tolist() for k, v in self.MetaObj.ViewingAngle_arrProv.items()},
meshwidth=10, # for faster processing self.shape_fullArr,
subset=None, meshwidth=10, # for faster processing
bandwise=0) subset=None,
bandwise=0)
else: else:
self._VZA_arr = GEOP.calc_VZA_array(self.shape_fullArr, self.arr_pos, self.fullSceneCornerPos, self._VZA_arr = GEOP.calc_VZA_array(self.shape_fullArr, self.arr_pos, self.fullSceneCornerPos,
float(self.meta_odict['ViewingAngle']), float(self.MetaObj.ViewingAngle),
float(self.meta_odict['FieldOfView']), float(self.MetaObj.FOV),
self.logger, self.logger,
nodata_mask=None, # dont overwrite areas outside image with nodata nodata_mask=None, # dont overwrite areas outside image with nodata
outFill=get_outFillZeroSaturated(np.float32)[0], outFill=get_outFillZeroSaturated(np.float32)[0],
...@@ -120,13 +121,14 @@ class L1C_object(L1B_object): ...@@ -120,13 +121,14 @@ class L1C_object(L1B_object):
""" """
if self._VAA_arr is None: if self._VAA_arr is None:
self.logger.info('Calculating viewing azimuth array...') self.logger.info('Calculating viewing azimuth array...')
if 'IncidenceAngle_arrProv' in self.meta_odict and self.meta_odict['IncidenceAngle_arrProv']: if self.MetaObj.IncidenceAngle_arrProv:
# Sentinel-2 # Sentinel-2
self._VAA_arr = GEOP.adjust_acquisArrProv_to_shapeFullArr(self.meta_odict['IncidenceAngle_arrProv'], self._VAA_arr = GEOP.adjust_acquisArrProv_to_shapeFullArr(
self.shape_fullArr, {k: v.tolist() for k, v in self.MetaObj.IncidenceAngle_arrProv.items()},
meshwidth=10, # for faster processing self.shape_fullArr,
subset=None, meshwidth=10, # for faster processing
bandwise=0) subset=None,
bandwise=0)
else: else:
# only a mean VAA is available # only a mean VAA is available
if self.VAA_mean is None: if self.VAA_mean is None:
...@@ -156,11 +158,11 @@ class L1C_object(L1B_object): ...@@ -156,11 +158,11 @@ class L1C_object(L1B_object):
self._SZA_arr, self._SAA_arr = \ self._SZA_arr, self._SAA_arr = \
GEOP.calc_SZA_SAA_array( GEOP.calc_SZA_SAA_array(
self.shape_fullArr, self.arr_pos, self.shape_fullArr, self.arr_pos,
self.meta_odict['AcqDate'], self.MetaObj.AcqDate,
self.meta_odict['AcqTime'], self.MetaObj.AcqTime,
self.fullSceneCornerPos, self.fullSceneCornerPos,
self.fullSceneCornerLonLat, self.fullSceneCornerLonLat,
self.meta_odict['overpass duraction sec'], self.MetaObj.overpassDurationSec,
self.logger, self.logger,
meshwidth=10, meshwidth=10,
nodata_mask=None, # dont overwrite areas outside the image with nodata nodata_mask=None, # dont overwrite areas outside the image with nodata
...@@ -294,7 +296,7 @@ class AtmCorr(object): ...@@ -294,7 +296,7 @@ class AtmCorr(object):
logger_atmCorr.addHandler(fileHandler) logger_atmCorr.addHandler(fileHandler)
inObj.close_GMS_loggers() inObj.close_loggers()
self._logger = logger_atmCorr self._logger = logger_atmCorr
return self._logger return self._logger
...@@ -315,7 +317,7 @@ class AtmCorr(object): ...@@ -315,7 +317,7 @@ class AtmCorr(object):
self._logger.close() self._logger.close()
self._logger = None self._logger = None
[inObj.close_GMS_loggers() for inObj in self.inObjs] [inObj.close_loggers() for inObj in self.inObjs]
@property @property
def GSDs(self): def GSDs(self):
...@@ -360,7 +362,7 @@ class AtmCorr(object): ...@@ -360,7 +362,7 @@ class AtmCorr(object):
# float32! -> conversion to np.float16 will convert -9999 to -10000 # float32! -> conversion to np.float16 will convert -9999 to -10000
arr2pass = inObj.arr[:, :, bandIdx].astype(np.float32) arr2pass = inObj.arr[:, :, bandIdx].astype(np.float32)
arr2pass[arr2pass == inObj.arr.nodata] = np.nan # set nodata values to np.nan arr2pass[arr2pass == inObj.arr.nodata] = np.nan # set nodata values to np.nan
data_dict[bandN] = (arr2pass / inObj.meta_odict['ScaleFactor']).astype(np.float16) data_dict[bandN] = (arr2pass / inObj.MetaObj.ScaleFactor).astype(np.float16)
else: else:
inObj.logger.warning("Band '%s' cannot be included into atmospheric correction because it " inObj.logger.warning("Band '%s' cannot be included into atmospheric correction because it "
"exists multiple times." % bandN) "exists multiple times." % bandN)
...@@ -472,14 +474,14 @@ class AtmCorr(object): ...@@ -472,14 +474,14 @@ class AtmCorr(object):
del self.logger # otherwise each input object would have multiple fileHandlers del self.logger # otherwise each input object would have multiple fileHandlers
metadata = dict( metadata = dict(
U=self.inObjs[0].meta_odict['EarthSunDist'], U=self.inObjs[0].MetaObj.EarthSunDist,
SENSING_TIME=self.inObjs[0].acq_datetime, SENSING_TIME=self.inObjs[0].acq_datetime,
# SENSING_TIME=datetime.strptime('2015-08-12 10:40:21 +0000', '%Y-%m-%d %H:%M:%S %z'), # SENSING_TIME=datetime.strptime('2015-08-12 10:40:21 +0000', '%Y-%m-%d %H:%M:%S %z'),
viewing_zenith=self._meta_get_viewing_zenith(), viewing_zenith=self._meta_get_viewing_zenith(),
viewing_azimuth=self._meta_get_viewing_azimuth(), viewing_azimuth=self._meta_get_viewing_azimuth(),
relative_viewing_azimuth=self._meta_get_relative_viewing_azimuth(), relative_viewing_azimuth=self._meta_get_relative_viewing_azimuth(),
sun_mean_azimuth=self.inObjs[0].meta_odict['SunAzimuth'], sun_mean_azimuth=self.inObjs[0].MetaObj.SunAzimuth,
sun_mean_zenith=90 - self.inObjs[0].meta_odict['SunElevation'], sun_mean_zenith=90 - self.inObjs[0].MetaObj.SunElevation,
solar_irradiance=self._meta_get_solar_irradiance(), solar_irradiance=self._meta_get_solar_irradiance(),
aux_data=self._meta_get_aux_data(), aux_data=self._meta_get_aux_data(),
spatial_samplings=self._meta_get_spatial_samplings() spatial_samplings=self._meta_get_spatial_samplings()
...@@ -676,8 +678,8 @@ class AtmCorr(object): ...@@ -676,8 +678,8 @@ class AtmCorr(object):
# FIXME calculation of center wavelengths within SRF() used not the GMS algorithm # FIXME calculation of center wavelengths within SRF() used not the GMS algorithm
# SRF instance must be created for all bands and the previous proc level # SRF instance must be created for all bands and the previous proc level
GMS_identifier_fullScene = self.inObjs[0].GMS_identifier GMS_identifier_fullScene = self.inObjs[0].GMS_identifier
GMS_identifier_fullScene['Subsystem'] = '' GMS_identifier_fullScene.subsystem = ''
GMS_identifier_fullScene['proc_level'] = proc_chain[proc_chain.index(self.inObjs[0].proc_level) - 1] GMS_identifier_fullScene.proc_level = proc_chain[proc_chain.index(self.inObjs[0].proc_level) - 1]
return SRF(GMS_identifier_fullScene, wvl_unit='nanometers', format_bandnames=True) return SRF(GMS_identifier_fullScene, wvl_unit='nanometers', format_bandnames=True)
...@@ -690,7 +692,7 @@ class AtmCorr(object): ...@@ -690,7 +692,7 @@ class AtmCorr(object):
tgt_res = self.inObjs[0].ac_options['cld_mask']['target_resolution'] tgt_res = self.inObjs[0].ac_options['cld_mask']['target_resolution']
# check if input GMS objects provide a cloud mask # check if input GMS objects provide a cloud mask
avail_cloud_masks = {inObj.GMS_identifier['Subsystem']: inObj.mask_clouds for inObj in self.inObjs} avail_cloud_masks = {inObj.GMS_identifier.subsystem: inObj.mask_clouds for inObj in self.inObjs}
no_avail_CMs = list(set(avail_cloud_masks.values())) == [None] no_avail_CMs = list(set(avail_cloud_masks.values())) == [None]
# compute cloud mask if not already provided # compute cloud mask if not already provided
...@@ -950,7 +952,6 @@ class AtmCorr(object): ...@@ -950,7 +952,6 @@ class AtmCorr(object):
inObj.MetaObj.LayerBandsAssignment = out_LBA inObj.MetaObj.LayerBandsAssignment = out_LBA
inObj.LayerBandsAssignment = out_LBA inObj.LayerBandsAssignment = out_LBA
inObj.MetaObj.filter_layerdependent_metadata() inObj.MetaObj.filter_layerdependent_metadata()
inObj.meta_odict = inObj.MetaObj.to_odict() # actually auto-updated by getter
################################################################################## ##################################################################################
# join SURFACE REFLECTANCE as 3D int16 array, scaled to scale factor from config # # join SURFACE REFLECTANCE as 3D int16 array, scaled to scale factor from config #
......
...@@ -36,6 +36,7 @@ from ..misc.definition_dicts import datasetid_to_sat_sen, sat_sen_to_datasetid ...@@ -36,6 +36,7 @@ from ..misc.definition_dicts import datasetid_to_sat_sen, sat_sen_to_datasetid
from ..misc.exceptions import ClassifierNotAvailableError from ..misc.exceptions import ClassifierNotAvailableError
from ..model.metadata import get_LayerBandsAssignment from ..model.metadata import get_LayerBandsAssignment
from .L2A_P import L2A_object from .L2A_P import L2A_object
from ..model.gms_object import GMS_identifier
__author__ = 'Daniel Scheffler' __author__ = 'Daniel Scheffler'
...@@ -60,13 +61,25 @@ class L2B_object(L2A_object): ...@@ -60,13 +61,25 @@ class L2B_object(L2A_object):
method = CFG.spechomo_method method = CFG.spechomo_method
src_dsID = sat_sen_to_datasetid(self.satellite, self.sensor) src_dsID = sat_sen_to_datasetid(self.satellite, self.sensor)
src_cwls = self.meta_odict['wavelength'] src_cwls = [float(self.MetaObj.CWL[bN]) for bN in self.MetaObj.LayerBandsAssignment]
# FIXME exclude or include thermal bands; respect sorted CWLs in context of LayerBandsAssignment # FIXME exclude or include thermal bands; respect sorted CWLs in context of LayerBandsAssignment
tgt_sat, tgt_sen = datasetid_to_sat_sen(CFG.datasetid_spectral_ref) tgt_sat, tgt_sen = datasetid_to_sat_sen(CFG.datasetid_spectral_ref)
# NOTE: get target LBA at L2A, because spectral characteristics of target sensor do not change after AC # NOTE: get target LBA at L2A, because spectral characteristics of target sensor do not change after AC
tgt_LBA = get_LayerBandsAssignment( tgt_LBA = get_LayerBandsAssignment(
dict(Satellite=tgt_sat, Sensor=tgt_sen, Subsystem=None, GMS_identifier(satellite=tgt_sat, sensor=tgt_sen, subsystem='',
image_type='RSD', proc_level='L2A', dataset_ID=src_dsID, logger=None)) image_type='RSD', proc_level='L2A', dataset_ID=src_dsID, logger=None))
if CFG.datasetid_spectral_ref is None:
tgt_cwl = CFG.target_CWL
tgt_fwhm = CFG.target_FWHM
else:
# exclude those bands from CFG.target_CWL and CFG.target_FWHM that have been removed after AC
full_LBA = get_LayerBandsAssignment(
GMS_identifier(satellite=tgt_sat, sensor=tgt_sen, subsystem='',
image_type='RSD', proc_level='L1A', dataset_ID=src_dsID, logger=None),
no_thermal=True, no_pan=False, return_fullLBA=True, sort_by_cwl=True, proc_level='L1A')
tgt_cwl = [dict(zip(full_LBA, CFG.target_CWL))[bN] for bN in tgt_LBA]
tgt_fwhm = [dict(zip(full_LBA, CFG.target_FWHM))[bN] for bN in tgt_LBA]
#################################################### ####################################################
# special cases where homogenization is not needed # # special cases where homogenization is not needed #
...@@ -92,7 +105,7 @@ class L2B_object(L2A_object): ...@@ -92,7 +105,7 @@ class L2B_object(L2A_object):
if method == 'LI' or CFG.datasetid_spectral_ref is None: 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) # linear interpolation (if intended by user or in case of custom spectral characteristics of target sensor)
# -> no classifier for that case available -> linear interpolation # -> no classifier for that case available -> linear interpolation
im = SpH.interpolate_cube(self.arr, src_cwls, CFG.target_CWL, kind='linear') im = SpH.interpolate_cube(self.arr, src_cwls, tgt_cwl, kind='linear')
if CFG.spechomo_estimate_accuracy: if CFG.spechomo_estimate_accuracy:
self.logger.warning("Unable to compute any error information in case spectral homogenization algorithm " self.logger.warning("Unable to compute any error information in case spectral homogenization algorithm "
...@@ -114,7 +127,7 @@ class L2B_object(L2A_object): ...@@ -114,7 +127,7 @@ class L2B_object(L2A_object):
compute_errors=CFG.spechomo_estimate_accuracy, compute_errors=CFG.spechomo_estimate_accuracy,
bandwise_errors=CFG.spechomo_bandwise_accuracy, bandwise_errors=CFG.spechomo_bandwise_accuracy,
fallback_argskwargs=dict( fallback_argskwargs=dict(
args=dict(source_CWLs=src_cwls, target_CWLs=CFG.target_CWL,), args=dict(source_CWLs=src_cwls, target_CWLs=tgt_cwl,),
kwargs=dict(kind='linear') kwargs=dict(kind='linear')
)) ))
...@@ -123,11 +136,9 @@ class L2B_object(L2A_object): ...@@ -123,11 +136,9 @@ class L2B_object(L2A_object):
################### ###################
self.LayerBandsAssignment = tgt_LBA self.LayerBandsAssignment = tgt_LBA
self.meta_odict['wavelength'] = list(CFG.target_CWL) self.MetaObj.CWL = dict(zip(tgt_LBA, tgt_cwl))
self.meta_odict['bandwidths'] = list(CFG.target_FWHM) self.MetaObj.FWHM = dict(zip(tgt_LBA, tgt_fwhm))
self.meta_odict['bands'] = len(CFG.target_CWL) self.MetaObj.bands = len(tgt_LBA)
if 'band names' in self.meta_odict: # FIXME bug workaround
del self.meta_odict['band names'] # TODO
self.arr = im # type: GeoArray self.arr = im # type: GeoArray
self.spec_homo_errors = errs # type: Union[np.ndarray, None] # int16, None if ~CFG.spechomo_estimate_accuracy self.spec_homo_errors = errs # type: Union[np.ndarray, None] # int16, None if ~CFG.spechomo_estimate_accuracy
...@@ -936,15 +947,15 @@ class ReferenceCube_Generator(object): ...@@ -936,15 +947,15 @@ class ReferenceCube_Generator(object):
return self._refcubes return self._refcubes
def _get_tgt_GMS_identifier(self, tgt_sat, tgt_sen): def _get_tgt_GMS_identifier(self, tgt_sat, tgt_sen):
# type: (str, str) -> dict # type: (str, str) -> GMS_identifier
"""Get a GMS identifier for the specified target sensor such that all possible bands are included (L1A) """Get a GMS identifier for the specified target sensor such that all possible bands are included (L1A)
:param tgt_sat: target satellite :param tgt_sat: target satellite
:param tgt_sen: target sensor :param tgt_sen: target sensor
:return: :return:
""" """
return dict(Satellite=tgt_sat, Sensor=tgt_sen, Subsystem=None, image_type='RSD', return GMS_identifier(satellite=tgt_sat, sensor=tgt_sen, subsystem=None, image_type='RSD', dataset_ID=None,
proc_level='L1A', logger=self.logger) # use L1A to have all bands available proc_level='L1A', logger=self.logger) # use L1A to have all bands available
def _get_tgt_LayerBandsAssignment(self, tgt_sat, tgt_sen): def _get_tgt_LayerBandsAssignment(self, tgt_sat, tgt_sen):
# type: (str, str) -> list # type: (str, str) -> list
...@@ -1035,8 +1046,8 @@ class ReferenceCube_Generator(object): ...@@ -1035,8 +1046,8 @@ class ReferenceCube_Generator(object):
# first, perform spectral resampling to Sentinel-2 to reduce dimensionality # first, perform spectral resampling to Sentinel-2 to reduce dimensionality
if downsamp_sat and downsamp_sen: if downsamp_sat and downsamp_sen:
tgt_srf = SRF(dict(Satellite=downsamp_sat, Sensor=downsamp_sen, Subsystem=None, image_type='RSD', tgt_srf = SRF(GMS_identifier(satellite=downsamp_sat, sensor=downsamp_sen, subsystem=None, image_type='RSD',
proc_level='L1A', logger=self.logger)) dataset_ID=None, proc_level='L1A', logger=self.logger))
im2clust = self.resample_image_spectrally(im2clust, tgt_srf, progress=progress) im2clust = self.resample_image_spectrally(im2clust, tgt_srf, progress=progress)