Commit a75ee7c3 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Added support for atmospheric correction of Landsat-7 and Landsat-8.

algorithms.gms_object.GMS_object:
- arr.setter: catched an error during setting of arr.bandnames
- revised ac_options

algorithms.L1C_P:
- L1C_object:
    - revised logger.deleter
    - added property 'options'
    - added _get_mask_clouds()
    - revised run_atmospheric_correction()
    - revised _join_results_to_inObjs()

database:
- removed folder ac_options -> AC options are now directly imported fom S2SCAPEM (thus version controlled)

misc.path_generator:
- revised get_path_ac_options() -> AC options are now directly imported fom S2SCAPEM (thus version controlled)

- updated __version__
parent 32d65014
...@@ -15,7 +15,7 @@ from . import config ...@@ -15,7 +15,7 @@ from . import config
from .processing.process_controller import process_controller from .processing.process_controller import process_controller
__version__ = '20170127.01' __version__ = '20170208.01'
__author__ = 'Daniel Scheffler' __author__ = 'Daniel Scheffler'
__all__ = ['algorithms', __all__ = ['algorithms',
'io', 'io',
......
...@@ -4,17 +4,7 @@ ...@@ -4,17 +4,7 @@
# Level 1C Processor: # Level 1C Processor:
# #
# Performed operations: # Performed operations:
# Atmospheric correction of radiance data: # Atmospheric correction of TOA-reflectance data:
# Inputs:
# - radiance array
# - DGM array
# - atmospheric layers (ozone,aerosols,water vapour)
# - sensor SRFs
# Outputs:
# - reflectance array
# - metadata dictionary with new information about:
# - quality of atmospheric correction
# - atmopsheric water vapour content, ...
# #
# Written by Daniel Scheffler # Written by Daniel Scheffler
# #
...@@ -43,7 +33,7 @@ from ..misc.definition_dicts import get_outFillZeroSaturated, proc_chain ...@@ -43,7 +33,7 @@ from ..misc.definition_dicts import get_outFillZeroSaturated, proc_chain
from ..io.Input_reader import SRF from ..io.Input_reader import SRF
from S2SCAPEM.s2_ac import AC_GMS from S2SCAPEM.s2_ac import AC_GMS
from S2MSI import RSImage from S2MSI import RSImage, S2Mask
########################### core functions #################################### ########################### core functions ####################################
...@@ -244,6 +234,7 @@ class AtmCorr(object): ...@@ -244,6 +234,7 @@ class AtmCorr(object):
self._metadata = {} self._metadata = {}
self._nodata = {} self._nodata = {}
self._band_spatial_sampling = {} self._band_spatial_sampling = {}
self._options = {}
# assertions # assertions
scene_IDs = [obj.scene_ID for obj in L1C_objs] scene_IDs = [obj.scene_ID for obj in L1C_objs]
...@@ -305,8 +296,9 @@ class AtmCorr(object): ...@@ -305,8 +296,9 @@ class AtmCorr(object):
@logger.deleter @logger.deleter
def logger(self): def logger(self):
self._logger.close() if self._logger not in [None, 'not set']:
self._logger = None self._logger.close()
self._logger = None
[inObj.close_GMS_loggers() for inObj in self.inObjs] [inObj.close_GMS_loggers() for inObj in self.inObjs]
...@@ -488,6 +480,16 @@ class AtmCorr(object): ...@@ -488,6 +480,16 @@ class AtmCorr(object):
return self._metadata return self._metadata
@property
def options(self):
"""Returns a dictionary containing AC options.
"""
if self._options:
return self._options
else:
self._options = self.inObjs[0].ac_options
return self._options
def _meta_get_spatial_samplings(self): def _meta_get_spatial_samplings(self):
""" """
...@@ -663,6 +665,46 @@ class AtmCorr(object): ...@@ -663,6 +665,46 @@ class AtmCorr(object):
return SRF(GMS_identifier_fullScene, wvl_unit='nanometers', format_bandnames=True) return SRF(GMS_identifier_fullScene, wvl_unit='nanometers', format_bandnames=True)
def _get_mask_clouds(self):
"""Returns an instance of S2Mask in case cloud mask is given by input GMS objects. Otherwise None is returned.
:return:
"""
# check if input GMS objects provide a cloud mask
avail_cloud_masks = {inObj.GMS_identifier['Subsystem']: inObj.mask_clouds for inObj in self.inObjs}
if list(set(avail_cloud_masks.values())) == [None]:
return None
# get cloud mask target resolution
tgt_res = self.inObjs[0].ac_options['cld_mask']['target_resolution']
inObjs2use = [obj for obj in self.inObjs if obj.mask_clouds is not None and obj.mask_clouds.xgsd==tgt_res]
if not inObjs2use:
raise ValueError('Error appending cloud mask to input arguments of atmospheric correction. No input '
'GMS object provides a cloud mask with spatial resolution of %s.' %tgt_res)
inObj2use = inObjs2use[0]
# get mask array
cm_array = inObj2use.masks[:,:,1] # FIXME hardcoded
# get legend
cm_legend = {'Clear': 10, 'Thick Clouds': 20, 'Thin Clouds': 30, 'Snow': 40} # FIXME hardcoded
# get geocoding
# validate GSD
if inObj2use.mask_clouds.xgsd != inObj2use.mask_clouds.ygsd:
warnings.warn("Cloud mask X/Y GSD is not equal for entity ID %s" % inObj2use.entity_ID +
(' (%s)' % inObj2use.subsystem if inObj2use.subsystem else '') +
'Using X-GSD as key for cloud mask geocoding.')
cm_geocoding = self.metadata["spatial_samplings"][tgt_res]
self.options['cld_mask']['nodata_value_mask'] = inObj2use.mask_clouds.nodata
return S2Mask(mask_array = cm_array,
mask_legend = cm_legend,
geo_coding = cm_geocoding)
def run_atmospheric_correction(self, dump_ac_input=False): def run_atmospheric_correction(self, dump_ac_input=False):
# type: (bool) -> list # type: (bool) -> list
"""Collects all input data for atmospheric correction, runs the AC and returns the corrected L1C objects """Collects all input data for atmospheric correction, runs the AC and returns the corrected L1C objects
...@@ -682,14 +724,18 @@ class AtmCorr(object): ...@@ -682,14 +724,18 @@ class AtmCorr(object):
band_spatial_sampling = self.band_spatial_sampling, band_spatial_sampling = self.band_spatial_sampling,
tile_name = self.tile_name, tile_name = self.tile_name,
dem = self._get_dem(), dem = self._get_dem(),
srf = self._get_srf() srf = self._get_srf(),
mask_clouds = self._get_mask_clouds() # returns an instance of S2Mask or None if cloud mask is not given by input GMS objects
) # NOTE: all keys of this dict are later converted to attributes of RSImage ) # NOTE: all keys of this dict are later converted to attributes of RSImage
options = self.inObjs[0].ac_options
script = False script = False
# create an instance of RSImage
rs_image = RSImage(**rs_data)
self.ac_input = dict( self.ac_input = dict(
rs_image = rs_data, rs_image = rs_image,
options = options, options = self.options,
logger = repr(self.logger), # only a string logger = repr(self.logger), # only a string
script = script script = script
) )
...@@ -697,9 +743,10 @@ class AtmCorr(object): ...@@ -697,9 +743,10 @@ class AtmCorr(object):
# run AC # run AC
self.logger.info('Atmospheric correction started.') self.logger.info('Atmospheric correction started.')
try: try:
self.results = AC_GMS(RSImage(**rs_data), options, logger=self.logger, script=False) rs_image.logger = self.logger
self.results = AC_GMS(rs_image, self.options, logger=self.logger, script=script)
except: except:
# serialialize AC inputs # serialialize AC input
if dump_ac_input: if dump_ac_input:
path_dump = self.inObjs[0].pathGen.get_path_ac_input_dump() path_dump = self.inObjs[0].pathGen.get_path_ac_input_dump()
with open(path_dump, 'wb') as outF: with open(path_dump, 'wb') as outF:
...@@ -715,7 +762,7 @@ class AtmCorr(object): ...@@ -715,7 +762,7 @@ class AtmCorr(object):
raise raise
# get processing infos # get processing infos
self.proc_info = self.ac_input['options']['processing'] self.proc_info = self.ac_input['options']['processing'] # FIXME this is not appended to GMS objects
# join results # join results
self._join_results_to_inObjs() # sets self.outObjs self._join_results_to_inObjs() # sets self.outObjs
...@@ -758,22 +805,33 @@ class AtmCorr(object): ...@@ -758,22 +805,33 @@ class AtmCorr(object):
inObj.arr = surf_refl.astype(inObj.arr.dtype) # -> int16 inObj.arr = surf_refl.astype(inObj.arr.dtype) # -> int16
# join ERRORS ARRAY as 3D int8 or int16 BOA reflectance array, scaled to scale factor from config # join ERRORS ARRAY as 3D int8 or int16 BOA reflectance array, scaled to scale factor from config
ac_errors = np.dstack((self.results.data_errors[bandN] for bandN in ac_bandNs)) if self.results.data_errors is not None:
ac_errors *= CFG.usecase.scale_factor_errors_ac # scale using scale factor (output is float16) ac_errors = np.dstack((self.results.data_errors[bandN] for bandN in ac_bandNs))
out_dtype = np.int8 if CFG.usecase.scale_factor_errors_ac<=255 else np.int16 ac_errors *= CFG.usecase.scale_factor_errors_ac # scale using scale factor (output is float16)
ac_errors[nodata] = get_outFillZeroSaturated(out_dtype)[0] out_dtype = np.int8 if CFG.usecase.scale_factor_errors_ac<=255 else np.int16
ac_errors = ac_errors.astype(out_dtype) ac_errors[nodata] = get_outFillZeroSaturated(out_dtype)[0]
inObj.ac_errors = ac_errors # setter generates a GeoArray with the same bandnames like inObj.arr ac_errors = ac_errors.astype(out_dtype)
# TODO how to handle nans? inObj.ac_errors = ac_errors # setter generates a GeoArray with the same bandnames like inObj.arr
# TODO how to handle nans?
else:
self.logger.warning("Atmospheric correction did not provide a 'data_errors' array. Maybe due to "
"missing SNR model? GMS_object.ac_errors kept None.")
# join CLOUD MASK as 2D uint8 array # join CLOUD MASK as 2D uint8 array
# NOTE: mask_clouds has also methods 'export_mask_rgb()', 'export_confidence_to_jpeg2000()', ... # NOTE: mask_clouds has also methods 'export_mask_rgb()', 'export_confidence_to_jpeg2000()', ...
mask_clouds_ac = self.results.mask_clouds.mask_array # uint8 2D array mask_clouds_ac = self.results.mask_clouds.mask_array # uint8 2D array
cfd_arr = self.results.mask_clouds.mask_confidence_array # float32 2D array, scaled [0-1, nodata 255]
cfd_arr[cfd_arr == self.ac_input['options']['cld_mask']['nodata_value_mask']] = -1 # join confidence array for mask_clouds
cfd_arr = (cfd_arr * CFG.usecase.scale_factor_BOARef).astype(np.int16) if self.results.mask_clouds.mask_confidence_array is not None:
cfd_arr[cfd_arr == -CFG.usecase.scale_factor_BOARef] = get_outFillZeroSaturated(cfd_arr.dtype)[0] cfd_arr = self.results.mask_clouds.mask_confidence_array # float32 2D array, scaled [0-1, nodata 255]
cfd_arr[cfd_arr == self.ac_input['options']['cld_mask']['nodata_value_mask']] = -1
cfd_arr = (cfd_arr * CFG.usecase.scale_factor_BOARef).astype(np.int16)
cfd_arr[cfd_arr == -CFG.usecase.scale_factor_BOARef] = get_outFillZeroSaturated(cfd_arr.dtype)[0]
else:
cfd_arr = None
self.logger.warning("Atmospheric correction did not provide a 'mask_confidence_array' array for "
"attribute 'mask_clouds. GMS_object.mask_clouds_confidence kept None.")
joined = False joined = False
for inObj in self.inObjs: for inObj in self.inObjs:
...@@ -800,8 +858,9 @@ class AtmCorr(object): ...@@ -800,8 +858,9 @@ class AtmCorr(object):
inObj.mask_clouds.nodata = mask_clouds_nodata inObj.mask_clouds.nodata = mask_clouds_nodata
# set cloud mask confidence array # set cloud mask confidence array
inObj.mask_clouds_confidence = GeoArray(cfd_arr, inObj.arr.gt, inObj.arr.prj, if cfd_arr is not None:
nodata=get_outFillZeroSaturated(cfd_arr.dtype)[0]) inObj.mask_clouds_confidence = GeoArray(cfd_arr, inObj.arr.gt, inObj.arr.prj,
nodata=get_outFillZeroSaturated(cfd_arr.dtype)[0])
joined=True joined=True
# update masks (always do that because masks can also only contain one layer) # update masks (always do that because masks can also only contain one layer)
...@@ -811,6 +870,9 @@ class AtmCorr(object): ...@@ -811,6 +870,9 @@ class AtmCorr(object):
self.logger.warning('Cloud mask has not been appended to one of the AC inputs because there was no' self.logger.warning('Cloud mask has not been appended to one of the AC inputs because there was no'
'input GMS object with the same dimensions.') 'input GMS object with the same dimensions.')
else:
self.logger.warning('Atmospheric correction did not return a result for the input array. '
'Thus the output keeps NOT atmospherically corrected.')
self.outObjs = self.inObjs self.outObjs = self.inObjs
......
...@@ -377,7 +377,11 @@ class GMS_object(object): ...@@ -377,7 +377,11 @@ class GMS_object(object):
# set bandnames like this: [B01, .., B8A,] # set bandnames like this: [B01, .., B8A,]
if self.LayerBandsAssignment: if self.LayerBandsAssignment:
self._arr.bandnames = self.LBA2bandnames(self.LayerBandsAssignment) if len(self.LayerBandsAssignment) == self._arr.bands:
self._arr.bandnames = self.LBA2bandnames(self.LayerBandsAssignment)
else:
warnings.warn("Cannot set 'bandnames' attribute of GMS_object.arr because LayerBandsAssignment has %s "
"bands and GMS_object.arr has %s bands." %(len(self.LayerBandsAssignment), self._arr.bands))
@arr.deleter @arr.deleter
...@@ -584,7 +588,7 @@ class GMS_object(object): ...@@ -584,7 +588,7 @@ class GMS_object(object):
if not self._ac_options: if not self._ac_options:
path_ac_options = PG.get_path_ac_options(self.GMS_identifier) path_ac_options = PG.get_path_ac_options(self.GMS_identifier)
if os.path.exists(path_ac_options): if path_ac_options and os.path.exists(path_ac_options):
opt_dict = get_ac_options(path_ac_options) opt_dict = get_ac_options(path_ac_options)
# update some file paths according to # update some file paths according to
...@@ -597,11 +601,15 @@ class GMS_object(object): ...@@ -597,11 +601,15 @@ class GMS_object(object):
opt_dict['cld_mask']['persistence_file'] = PG.get_path_cloud_class_obj(self.GMS_identifier) opt_dict['cld_mask']['persistence_file'] = PG.get_path_cloud_class_obj(self.GMS_identifier)
opt_dict['output'] = [] # outputs are not needed for GMS -> so opt_dict['output'] = [] # outputs are not needed for GMS -> so
opt_dict['report']['report_path'] = os.path.join(self.pathGen.get_path_procdata(), '[TYPE]') opt_dict['report']['report_path'] = os.path.join(self.pathGen.get_path_procdata(), '[TYPE]')
opt_dict['uncertainties']['snr_model'] = PG.get_path_snr_model(self.GMS_identifier) if 'uncertainties' in opt_dict:
opt_dict['uncertainties']['snr_model'] = PG.get_path_snr_model(self.GMS_identifier)
#opt_dict['AC']['n_cores'] = CFG.job.CPUs if CFG.job.allow_subMultiprocessing else 1 #opt_dict['AC']['n_cores'] = CFG.job.CPUs if CFG.job.allow_subMultiprocessing else 1
self._ac_options = opt_dict self._ac_options = opt_dict
else:
self.logger.warning('There is no options file available for atmospheric correction. '
'Atmospheric correction must be skipped.')
return self._ac_options return self._ac_options
......
{
"ram": {
"upper_limit": 20,
"unit": "GB"
},
"ozo_to_DU": 71524.3,
"ozone_fallback":500.0, /* in [DU] */
"dT/dh":-0.006, /* in [K/m]*/
"S2Image": {
"S2_MSI_granule_path": "None",
"sliceY": "slice(None, None, None)",
"sliceX": "slice(None, None, None)",
"unit": "reflectance",
"aux_fields": {
"spr": [
20.0
],
"cwv": [
60.0,
20.0
],
"ozo": "mean"
},
"target_resolution": "None",
"import_bands": "all",
"driver": "gdal_JP2KAK" /*gdal_JP2KAK or gdal_JP2OpenJPEG*/
},
"ECMWF": {
"total_AOT": "fc_total_AOT_550nm",
"path_db": "/home/gts2/data/aux_data/ECMWF/",
"variables": [
"fc_total_AOT_550nm",
"fc_sulphate_AOT_550nm",
"fc_black_carbon_AOT_550nm",
"fc_dust_AOT_550nm",
"fc_organic_matter_AOT_550nm",
"fc_sea_salt_AOT_550nm"
],
"var2type": {
"fc_organic_matter_AOT_550nm": "aerosol_2",
"fc_sulphate_AOT_550nm": "aerosol_2",
"fc_sea_salt_AOT_550nm": "aerosol_2",
"fc_black_carbon_AOT_550nm": "aerosol_2",
"fc_dust_AOT_550nm": "aerosol_2"
}
},
"base_output_types": ["L2A","rgb_jpeg","metadata"],
"output": [],
"logger": {
"pprint": {
"indent": 3,
"compact": false,
"width": 100
},
"logg_options": true,
"datefmt": "%Y%m%d-%H:%M:%S",
"level": "INFO",
"format": "%(asctime)s:%(levelname)s::%(message)s",
"name": "S2SCAPEM"
},
"DEM": {
"fn": "/home/gts2/data/aux_data/dem/GLOBAL_SRTM_90m.h5",
"return_zeros_if_missing": true
},
"RTFO": {
"aerosol_1": {
"flag": 10,
"table_path": "/table_aerosol/type_1",
"dim_scat": [
"tau_a"
],
"atm_tables_fn": "/usr/local/gts2/linear_atm_functions_ncwv_5_npre_4_ncoz_2_ntmp_2_wvl_350.0_2550.0_1.00_pca.h5",
"dim_atm": [
"spr",
"coz",
"cwv",
"tmp"
],
"only_toa": true
},
"aerosol_0": {
"flag": 10,
"table_path": "/table_aerosol/type_0",
"dim_scat": [
"tau_a"
],
"atm_tables_fn": "/usr/local/gts2/linear_atm_functions_ncwv_5_npre_4_ncoz_2_ntmp_2_wvl_350.0_2550.0_1.00_pca.h5",
"dim_atm": [
"spr",
"coz",
"cwv",
"tmp"
],
"only_toa": true
},
"cirrus": {
"flag": 20,
"table_path": "/table_cirrus",
"dim_scat": [
"tau_c",
"eff_rad_cirris"
],
"atm_tables_fn": "/usr/local/gts2/linear_atm_functions_ncwv_5_npre_4_ncoz_2_ntmp_2_wvl_350.0_2550.0_1.00_pca.h5",
"dim_atm": [
"spr",
"coz",
"cwv",
"tmp"
],
"only_toa": true
},
"aerosol_2": {
"flag": 10,
"table_path": "/table_aerosol/type_2",
"dim_scat": [
"tau_a"
],
"atm_tables_fn": "/usr/local/gts2/linear_atm_functions_ncwv_5_npre_4_ncoz_2_ntmp_2_wvl_350.0_2550.0_1.00_pca.h5",
"dim_atm": [
"spr",
"coz",
"cwv",
"tmp"
],
"only_toa": true
},
"aerosol_3": {
"flag": 10,
"table_path": "/table_aerosol/type_3",
"dim_scat": [
"tau_a"
],
"atm_tables_fn": "/usr/local/gts2/linear_atm_functions_ncwv_5_npre_4_ncoz_2_ntmp_2_wvl_350.0_2550.0_1.00_pca.h5",
"dim_atm": [
"spr",
"coz",
"cwv",
"tmp"
],
"only_toa": true
}
},
"report": {
"reporting": false, /* either true or false */
"RGB": {
"output_size": 5490,
"hist_chop_off_fraction": 0.05
},
"HTML": true,
"figs": {
"MSK_rgb": {},
"DEM": {},
"CWV": {},
"AOT": {},
"RGB": {},
"SPR": {}
},
"figs_clouds": {
"MSK_rgb": {},
"DEM": {},
"RGB": {}
},
"report_path": "/home/gts2/logs/L2A/processing_reports/[TYPE]/",
"JPG": true,
"dpi": 200,
"n_cols": 2
},
"cld_mask": {
"target_resolution": 20.0,
"nodata_value_mask":255,
"persistence_file": "/home/gts2/data/aux_data/cld_mask/cld_mask_20160321_s2_8l0.200_11l0.060_12l0.040_v20161124_16:41:04.h5", /*base version: /home/AC/cld_mask_20160321_s2.h5 */
"processing_tiles": 4,
"majority_mask_filter": {
"majority_width": 5,
"block_replace": true,
"block_replace_params":[[50, 5]]
}
},
"AC": {
"aerosol_type_estimation":"maximum_ECMWF_type", /*should be:maximum_ECMWF_type*/
"max_pixel_processing": 30140100,
"n_smpl": {
"spr": 2,
"cwv": 3,
"tau_a": 2,
"rho": 10
},
"min_clear_fraction": 0.1,
"override_aerosol_type": "auto",
"parameter_bounds": "image",
"n_cores": 1,
"bands": [
"B01",
"B02",
"B03",
"B04",
"B05",
"B06",
"B07",
"B08",
"B8A",
"B11",
"B12"
]
},
"water_vapor": {
"int_order": 3,
"nanmean": 20,
"n_smpl": {
"spr": 5,
"cwv": 10,
"tau_a": 5,
"rho": 20
},
"sigma": 20,
"n_std": 1,
"bands": [
"B8A",
"B09"
],
"target_resolution": 60
},
"SRF": "S2A",
"uncertainties":{"snr_model":"/home/gts2/data/aux_data/S2A/S2A_SNR_model.csv"},
"aot_retrieval": {
"n_smpl": {
"spr": 2,
"tau_a": 5,
"rho": 5
},
"dark_nbins": 200,
"dark_fraction": "0.01",
"dark_bands": [
"B01",
"B02"
],
"dark_n_target": 500,
"bands": [
"B01",
"B02",
"B03",
"B04"
]
}
}
\ No newline at end of file
...@@ -252,7 +252,7 @@ class SRF(object): ...@@ -252,7 +252,7 @@ class SRF(object):
def from_dict(self, srf_dict): def from_dict(self, srf_dict):
"""Create an instance of SRF from a dictionary. """Create an instance of SRF from a dictionary.
:param dict: {'key_LayerBandsAssignment': <2D array: cols=[wvl,resp],rows=samples>} :param srf_dict: {'key_LayerBandsAssignment': <2D array: cols=[wvl,resp],rows=samples>}
""" """
is_nm = [300 <