Commit 0657813b 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__


Former-commit-id: a75ee7c3
parent dc5b1c37
......@@ -15,7 +15,7 @@ from . import config
from .processing.process_controller import process_controller
__version__ = '20170127.01'
__version__ = '20170208.01'
__author__ = 'Daniel Scheffler'
__all__ = ['algorithms',
'io',
......
......@@ -4,17 +4,7 @@
# Level 1C Processor:
#
# Performed operations:
# Atmospheric correction of radiance 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, ...
# Atmospheric correction of TOA-reflectance data:
#
# Written by Daniel Scheffler
#
......@@ -43,7 +33,7 @@ from ..misc.definition_dicts import get_outFillZeroSaturated, proc_chain
from ..io.Input_reader import SRF
from S2SCAPEM.s2_ac import AC_GMS
from S2MSI import RSImage
from S2MSI import RSImage, S2Mask
########################### core functions ####################################
......@@ -244,6 +234,7 @@ class AtmCorr(object):
self._metadata = {}
self._nodata = {}
self._band_spatial_sampling = {}
self._options = {}
# assertions
scene_IDs = [obj.scene_ID for obj in L1C_objs]
......@@ -305,8 +296,9 @@ class AtmCorr(object):
@logger.deleter
def logger(self):
self._logger.close()
self._logger = None
if self._logger not in [None, 'not set']:
self._logger.close()
self._logger = None
[inObj.close_GMS_loggers() for inObj in self.inObjs]
......@@ -488,6 +480,16 @@ class AtmCorr(object):
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):
"""
......@@ -663,6 +665,46 @@ class AtmCorr(object):
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):
# type: (bool) -> list
"""Collects all input data for atmospheric correction, runs the AC and returns the corrected L1C objects
......@@ -682,14 +724,18 @@ class AtmCorr(object):
band_spatial_sampling = self.band_spatial_sampling,
tile_name = self.tile_name,
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
options = self.inObjs[0].ac_options
script = False
# create an instance of RSImage
rs_image = RSImage(**rs_data)
self.ac_input = dict(
rs_image = rs_data,
options = options,
rs_image = rs_image,
options = self.options,
logger = repr(self.logger), # only a string
script = script
)
......@@ -697,9 +743,10 @@ class AtmCorr(object):
# run AC
self.logger.info('Atmospheric correction started.')
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:
# serialialize AC inputs
# serialialize AC input
if dump_ac_input:
path_dump = self.inObjs[0].pathGen.get_path_ac_input_dump()
with open(path_dump, 'wb') as outF:
......@@ -715,7 +762,7 @@ class AtmCorr(object):
raise
# 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
self._join_results_to_inObjs() # sets self.outObjs
......@@ -758,22 +805,33 @@ class AtmCorr(object):
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
ac_errors = np.dstack((self.results.data_errors[bandN] for bandN in ac_bandNs))
ac_errors *= CFG.usecase.scale_factor_errors_ac # scale using scale factor (output is float16)
out_dtype = np.int8 if CFG.usecase.scale_factor_errors_ac<=255 else np.int16
ac_errors[nodata] = get_outFillZeroSaturated(out_dtype)[0]
ac_errors = ac_errors.astype(out_dtype)
inObj.ac_errors = ac_errors # setter generates a GeoArray with the same bandnames like inObj.arr
# TODO how to handle nans?
if self.results.data_errors is not None:
ac_errors = np.dstack((self.results.data_errors[bandN] for bandN in ac_bandNs))
ac_errors *= CFG.usecase.scale_factor_errors_ac # scale using scale factor (output is float16)
out_dtype = np.int8 if CFG.usecase.scale_factor_errors_ac<=255 else np.int16
ac_errors[nodata] = get_outFillZeroSaturated(out_dtype)[0]
ac_errors = ac_errors.astype(out_dtype)
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
# 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
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]
# join confidence array for mask_clouds
if self.results.mask_clouds.mask_confidence_array is not None:
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
for inObj in self.inObjs:
......@@ -800,8 +858,9 @@ class AtmCorr(object):
inObj.mask_clouds.nodata = mask_clouds_nodata
# set cloud mask confidence array
inObj.mask_clouds_confidence = GeoArray(cfd_arr, inObj.arr.gt, inObj.arr.prj,
nodata=get_outFillZeroSaturated(cfd_arr.dtype)[0])
if cfd_arr is not None:
inObj.mask_clouds_confidence = GeoArray(cfd_arr, inObj.arr.gt, inObj.arr.prj,
nodata=get_outFillZeroSaturated(cfd_arr.dtype)[0])
joined=True
# update masks (always do that because masks can also only contain one layer)
......@@ -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'
'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
......
......@@ -377,7 +377,11 @@ class GMS_object(object):
# set bandnames like this: [B01, .., B8A,]
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
......@@ -584,7 +588,7 @@ class GMS_object(object):
if not self._ac_options:
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)
# update some file paths according to
......@@ -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['output'] = [] # outputs are not needed for GMS -> so
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
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
......
{
"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):
def from_dict(self, srf_dict):
"""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 < np.max(srf_dict[band][:, 0]) < 15000 for band in srf_dict]
......
......@@ -176,7 +176,7 @@ def get_path_cloud_class_obj(GMS_identifier, get_all=False):
"""
GMS_sensorcode = get_GMS_sensorcode(GMS_identifier)
satellite,sensor,logger = (GMS_identifier['Satellite'],GMS_identifier['Sensor'],GMS_identifier['logger'])
satellite,sensor,logger = GMS_identifier['Satellite'],GMS_identifier['Sensor'],GMS_identifier['logger']
path_cloud_classifier_objects = os.path.join(CFG.job.path_cloud_classif,satellite,sensor)
obj_name_dic = {
......@@ -279,14 +279,49 @@ def get_path_snr_model(GMS_identifier):
def get_path_ac_options(GMS_identifier):
# type: (dict)->str
# type: (dict)->any
"""Returns the path of the options json file needed for atmospheric correction.
"""
satellite,sensor = (GMS_identifier['Satellite'],GMS_identifier['Sensor'])
satellite = 'RapidEye' if re.match('RapidEye',satellite,re.I) else satellite
sensor = sensor[:-1] if re.match('SPOT',satellite,re.I) and sensor[-1] not in ['1','2'] else sensor
return os.path.join(CFG.job.path_ac_options, satellite, sensor, 'ac_options.json')
GMSid_ac = GMS_identifier
GMSid_ac['Subsystem'] = ''
sensorcode = get_GMS_sensorcode(GMSid_ac)
ac_options_file_dic = {
'AVNIR-2': None,
'TM4': None, # 'l8_options.json' # not yet usable due to missing ECMWF archive < 2012
'TM5': None, # 'l8_options.json' # not yet usable due to missing ECMWF archive < 2012
'TM7': 'l8_options.json', # AC uses Landsat-8 options for L7 but reads only a subset of the options
'LDCM': 'l8_options.json',
'SPOT1a': None,
'SPOT1b': None,
'SPOT2a': None,
'SPOT2b': None,
'SPOT3a': None,
'SPOT3b': None,
'SPOT4a': None,
'SPOT4b': None,
'SPOT5a': None,
'SPOT5b': None,
'RE5': None,
'AST_full': None,
'S2A_full': 's2_options.json',
'S2B_full': 's2_options.json',
}
try:
fName_optFile = ac_options_file_dic[get_GMS_sensorcode(GMS_identifier)]
except KeyError:
GMS_identifier['logger'].warning(
"Sensorcode '%s' is not included in ac_options dictionary. "
"Thus atmospheric correction is not available for the current scene." %sensorcode)
fName_optFile = None
if fName_optFile:
import S2SCAPEM
return os.path.join(os.path.dirname(S2SCAPEM.__file__), fName_optFile)
else:
return None
def get_path_ac_table(ac_options_key):
......
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