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

second (not completely working) version of wrapper for atmospheric correction

algorithms.GEOPROCESSING:
- revised imports
algorithms.gms_object:
- gms_object:
    - revised imports
    - added property 'dem': gms_object can now directly provide a corresponding SRTM DEM
    - arr: bandnames of property 'arr' are now in the form [B01, .., B8A,] and correspond to LayerBandsAssignment
    - added property 'ac_options': getter or options dictionary needed for atmospheric correction
    - from_disk(): added time zone to 'acquisition_date' datetime object
- added class failed_GMS_object (based on earlier version from helper functions)
algorithms.L1A_P.L1A_object:
- revised imports
- calc_TOARadRefTemp(): bugfix for wrong nodata value in out returned array in case of Sentinel-2
- update_spec_vals_according_to_dtype: bugfix for not updating L1A_object.arr.nodata
algorithms.L1B_P:
- revised imports
algorithms.L1C_P:
- L1C_object:
    - get_lonlat_coord_array(): changed handling of return values
    - calc_acquisition_illumination_geometry(): changed handling of return values
- AtmCorr:
    -  added attribute 'ac_input' containing input args/kwargs of atmospheric correction
    - data: now uses L1C_obj.arr.bandnames for lopping over bands
    - added property 'nodata'
    - added property 'tile_name'
    - added property 'band_spatial_sampling'
    - added property 'nodata'
    - added property 'nodata'
    - revised property 'metadata'
    - added _meta_get_spatial_samplings()
    - added _meta_get_solar_irradiance()
    - added _meta_get_viewing_zenith()
    - added _meta_get_viewing_azimuth()
    - added _meta_get_relative_viewing_azimuth()
    - revised run_atmospheric_correction()
    - revised join_results_to_inObjs()
algorithms.METADATA
- revised imports
- added 'ScaleFactor' to meta_odict
io.Input_reader:
- fixed some bad type hints
- SRF_reader(): moved path generator functionality to path_generator
- added open_specific_file_within_archive() (moved)
- added get_dem_by_extent(): new function for reading SRTM DEM data and warping to a given pixel grid
io.Output_writer:
- revised imports
- fixed some bad type hints
- added 'ScaleFactor' to enviHdr_keyOrder
misc.__init__:
- added __all__
misc.database_tools:
- fixed some bad type hints
- get_overlapping_scenes_from_postgreSQLdb(): bugfix for wrong indexing
misc.definition_dicts:
- new module, consisting of earlier functions from helper_functions
misc.exception_handler:
- new module, consisting of earlier functions from helper_functions
misc.helper_functions:
- moved trace_unhandled_exceptions(), log_uncaught_exceptions() to misc.exception_handler
- moved failed_GMS_object to gms_object
- moved get_job_summary to process_controller
- fixed some bad type hints
- moved get_GMS_sensorcode(), get_mask_classdefinition(), get_outFillZeroSaturated(), get_mask_colormap() to misc.definition_dicts
- moved open_specific_file_within_archive() to Input_reader
misc.path_generator:
- path_generator:
    - revised get_path_rawdata()
    - revised get_local_archive_path_baseN(): 'image_type' is not used anymore; removed deprecated warning
    - get_path_cloud_class_obj(): added cloud classificator files for Sentinel-2
    - added get_path_srf_file()
    - added get_path_snr_model()
    - added get_path_ac_options()
    - added get_path_ac_table()
processing.multiproc:
- revised MAP(): added new keyword 'flatten_ouput'
processing.pipeline:
- updated imports
- revised L1C_map(): input represents one OR multiple L1B_objects belonging to the same scene_ID (atmospheric correction has to be applied to ALL subsystems of a scene at once)
processing.process_controller:
- changed some map calls due to new keyword  'flatten_ouput' of processing.multiproc.MAP
- revised L1C_processing(): added grouping of L1B_objects by scene_ID
- revised create_job_summary(): bugfix for emtpy input list of get_job_summary()
- added get_job_summary(): moved from helper_functions
config.Job:
- added new attributes 'path_ac_options', 'path_ac_tables', 'path_SNR_models', 'path_dem_proc_srtm_90m', 'path_ECMWF_db'
pgSGL_db table 'config':
- added keys 'path_SNR_models', 'path_ac_options', 'path_dem_proc_srtm_90m', 'path_ECMWF_db',  'path_ac_tables'
Fileserver:
- added input datasets for atmospheric correction
- updated __version__
parent dd8f8334
......@@ -15,7 +15,7 @@ from . import config
from .processing.process_controller import process_controller
__version__ = '20170106.01'
__version__ = '20170111.01'
__author__ = 'Daniel Scheffler'
__all__ = ['algorithms',
'io',
......
......@@ -56,8 +56,10 @@ from py_tools_ds.ptds.geo.projection import get_UTMzone, EPSG2WKT, isProje
from py_tools_ds.ptds.geo.raster.reproject import warp_ndarray
from ..config import GMS_config as CFG
from ..io import envifilehandling as ef
from ..misc import helper_functions as HLP_F
from ..io import envifilehandling as ef
from ..misc.definition_dicts import get_outFillZeroSaturated
class GEOPROCESSING(object):
......@@ -390,7 +392,7 @@ class GEOPROCESSING(object):
self.logger.info(("\t*outPath: %s\n" % (self.workspace if outPath is None or outPath == "" else outPath)))
# --Define default for Special Values of reflectanceband
outFill, outZero, outSaturated = HLP_F.get_outFillZeroSaturated('int16')
outFill, outZero, outSaturated = get_outFillZeroSaturated('int16')
# --Create Outputfile
# Define name of outputfile
......@@ -986,8 +988,8 @@ class GEOPROCESSING(object):
t0= time.time()
in_arr = np.swapaxes(np.swapaxes(gdalnumeric.LoadFile(self.desc),0,2),0,1)
print('reading time', time.time()-t0)
if inFill is None: inFill = HLP_F.get_outFillZeroSaturated(in_arr.dtype)[0]
out_nodataVal = HLP_F.get_outFillZeroSaturated(in_arr.dtype)[0]
if inFill is None: inFill = get_outFillZeroSaturated(in_arr.dtype)[0]
out_nodataVal = get_outFillZeroSaturated(in_arr.dtype)[0]
t0= time.time()
warped, gt, prj = warp_ndarray(in_arr, self.geotransform, self.projection, out_prj,
......@@ -1010,8 +1012,8 @@ class GEOPROCESSING(object):
"""needs temporary files but does support multiprocessing and configuring cache size"""
t0 = time.time()
in_dtype = gdalnumeric.LoadFile(self.desc,0,0,1,1).dtype
if inFill is None: inFill = HLP_F.get_outFillZeroSaturated(in_dtype)[0]
out_nodataVal = HLP_F.get_outFillZeroSaturated(in_dtype)[0]
if inFill is None: inFill = get_outFillZeroSaturated(in_dtype)[0]
out_nodataVal = get_outFillZeroSaturated(in_dtype)[0]
gcps = ' '.join(['-gcp %s %s %s %s' %tuple(gcp) for gcp in [gcp_ul,gcp_ur,gcp_ll,gcp_lr]])
if use_workspace:
......@@ -1103,7 +1105,7 @@ class GEOPROCESSING(object):
:param custom_nodataVal:
"""
nodataVal = HLP_F.get_outFillZeroSaturated(np.int16)[0] if custom_nodataVal is None else custom_nodataVal
nodataVal = get_outFillZeroSaturated(np.int16)[0] if custom_nodataVal is None else custom_nodataVal
in_arr = array if array is not None else np.swapaxes(np.swapaxes(gdalnumeric.LoadFile(self.desc),0,2),0,1)
return np.all(np.where(in_arr==nodataVal,0,1),axis=2)
......@@ -1115,7 +1117,7 @@ class GEOPROCESSING(object):
mask_1bit = np.zeros((self.rows, self.cols, self.bands), dtype='uint8')
# img = np.empty((self.rows, self.cols, self.bands), dtype='int16')
nodataVal = HLP_F.get_outFillZeroSaturated(np.int16)[0] if custom_nodataVal is None else custom_nodataVal
nodataVal = get_outFillZeroSaturated(np.int16)[0] if custom_nodataVal is None else custom_nodataVal
# self.bandsList == range(self.bands) if subset != custom!
for x, out_idx in zip(self.bandsList, range(self.bands)):
band = self.inDs.GetRasterBand(x + 1) # Read inputband as numpyarray
......@@ -2813,9 +2815,9 @@ def DN2Rad(ndarray, offsets, gains, inFill=None, inZero=None, inSaturated=None,
"DN2Rad: Expected float or integer values from argument '%s'. Got %s" % (argname,type(i))
if [inFill,inZero,inSaturated] == [None,None,None]:
inFill,inZero,inSaturated = HLP_F.get_outFillZeroSaturated(ndarray.dtype)
inFill,inZero,inSaturated = get_outFillZeroSaturated(ndarray.dtype)
# --Define default for Special Values of reflectanceband
outFill, outZero, outSaturated = HLP_F.get_outFillZeroSaturated('int16')
outFill, outZero, outSaturated = get_outFillZeroSaturated('int16')
off,gai = [np.array(param).reshape(1,1,bands) for param in [offsets,gains]]
rad = (10*(off+ndarray*gai)).astype(np.int16)
......@@ -2863,9 +2865,9 @@ def DN2TOARef(ndarray, offsets, gains, irradiances, zenith, earthSunDist,
"DN2TOARef: Expected float or integer values from argument '%s'. Got %s" % (argname,type(i))
if [inFill,inZero,inSaturated] == [None,None,None]:
inFill,inZero,inSaturated = HLP_F.get_outFillZeroSaturated(ndarray.dtype)
inFill,inZero,inSaturated = get_outFillZeroSaturated(ndarray.dtype)
# --Define default for Special Values of reflectanceband
outFill, outZero, outSaturated = HLP_F.get_outFillZeroSaturated('int16')
outFill, outZero, outSaturated = get_outFillZeroSaturated('int16')
constant = 10000*math.pi*earthSunDist**2/math.cos(math.radians(zenith))
off,gai,irr = [np.array(param).reshape(1,1,bands) for param in [offsets,gains,irradiances]]
......@@ -2904,9 +2906,9 @@ def TOARad2Kelvin_fastforward(ndarray, K1, K2, emissivity=0.95, inFill=None, inZ
if [inFill,inZero,inSaturated] == [None,None,None]:
inFill,inZero,inSaturated = HLP_F.get_outFillZeroSaturated(ndarray.dtype)
inFill,inZero,inSaturated = get_outFillZeroSaturated(ndarray.dtype)
# --Define default for Special Values of reflectanceband
outFill, outZero, outSaturated = HLP_F.get_outFillZeroSaturated('int16')
outFill, outZero, outSaturated = get_outFillZeroSaturated('int16')
K1,K2 = [np.array(param).reshape(1,1,bands) for param in [K1,K2]]
Kelvin = (K2 / np.log(K1 / ndarray + 1)).astype(np.int16)
......@@ -2946,9 +2948,9 @@ def DN2DegreesCelsius_fastforward(ndarray,offsets, gains, K1, K2, emissivity=0.9
values from argument emissivity. Got %s" % type(emissivity)
if [inFill,inZero,inSaturated] == [None,None,None]:
inFill,inZero,inSaturated = HLP_F.get_outFillZeroSaturated(ndarray.dtype)
inFill,inZero,inSaturated = get_outFillZeroSaturated(ndarray.dtype)
# --Define default for Special Values of reflectanceband
outFill, outZero, outSaturated = HLP_F.get_outFillZeroSaturated('int16')
outFill, outZero, outSaturated = get_outFillZeroSaturated('int16')
off,gai,K1,K2 = [np.array(param).reshape(1,1,bands) for param in [offsets,gains,K1,K2]]
degCel = (100*((K2/np.log(K1*emissivity/(off+ndarray*gai)+1))-273.15)).astype(np.int16)
......@@ -3017,7 +3019,7 @@ def get_lonlat_coord_array(shape_fullArr, arr_pos, geotransform, projection, nod
xcoord_ycoord_arr = np.dstack((xcoord_arr, ycoord_arr))
if nodata_mask is not None:
outFill = outFill if outFill else HLP_F.get_outFillZeroSaturated('float32')[0]
outFill = outFill if outFill else get_outFillZeroSaturated('float32')[0]
xcoord_ycoord_arr[nodata_mask.astype(np.int8)==False] = outFill
UL_lonlat = (round(xcoord_ycoord_arr[0, 0, 0], 10), round(xcoord_ycoord_arr[0, 0, 1], 10))
......@@ -3084,7 +3086,7 @@ def calc_VZA_array(shape_fullArr, arr_pos, trueDataCornerPos, viewing_angle, FOV
cols_arr * rows_arr).astype(np.float32)
if nodata_mask is not None:
outFill = outFill if outFill else HLP_F.get_outFillZeroSaturated('float32')[0]
outFill = outFill if outFill else get_outFillZeroSaturated('float32')[0]
VZA_array[nodata_mask.astype(np.int8) == False] = outFill
return VZA_array
......
......@@ -27,6 +27,7 @@ from . import gms_cloud_classifier as CLD_P # Cloud Processor
from ..io import Output_writer as OUT_W
from ..misc import helper_functions as HLP_F
from ..misc import path_generator as PG
from ..misc.definition_dicts import get_outFillZeroSaturated
from .gms_object import GMS_object
......@@ -549,14 +550,20 @@ class L1A_object(GMS_object):
res = inArray
self.log_for_fullArr_or_firstTile('The input data already represents TOA '
'reflectance with the aimed scale factor of %d.' %CFG.usecase.scale_factor_TOARef)
else: # arr_desc == 'Temp'
raise NotImplementedError("Conversion Temp to %s is currently not supported." %conv)
# ensure int16 as output data type (also relevant for auto-setting of nodata value)
if isinstance(res,(np.ndarray, GeoArray)):
# change data type to int16 and update nodata values within array
res = res if res.dtype==np.int16 else res.astype(np.int16)
res[res==inFill] = get_outFillZeroSaturated(np.int16)[0]
if isinstance(res,np.ndarray): res = res if res.dtype==np.int16 else res.astype(np.int16)
if optical_thermal == 'optical': data_optical, optical_bandsList = res, bList
else: data_thermal, thermal_bandsList = res, bList
if optical_thermal == 'optical':
data_optical, optical_bandsList = res, bList
else:
data_thermal, thermal_bandsList = res, bList
# combine results from optical and thermal data
if data_optical is not None and data_thermal is not None:
......@@ -571,7 +578,7 @@ class L1A_object(GMS_object):
else:
dataOut = data_optical if data_thermal is None else data_thermal
assert dataOut is not None
#
self.update_spec_vals_according_to_dtype('int16')
tiles_desc = '_'.join([desc for op_th,desc in zip(['optical','thermal'], [CFG.usecase.conversion_type_optical,
CFG.usecase.conversion_type_thermal]) if desc in self.dict_LayerOptTherm.values()])
......@@ -756,7 +763,7 @@ class L1A_object(GMS_object):
self.mask_clouds = mask_clouds.astype(np.uint8)
# update fill values (fill values written by CLD_obj(self) are not equal to GMS fill values)
outFill = HLP_F.get_outFillZeroSaturated(self.mask_clouds.dtype)[0]
outFill = get_outFillZeroSaturated(self.mask_clouds.dtype)[0]
if outFill in mask_clouds:
warnings.warn('Value %s is assigned to mask_clouds although it already contains this value. '
......@@ -853,8 +860,10 @@ class L1A_object(GMS_object):
gdalnumeric.LoadFile(self.MetaObj.Dataname,0,0,1,1)
assert arr is not None
dtype = arr.dtype
self.MetaObj.spec_vals['fill'],self.MetaObj.spec_vals['zero'],self.MetaObj.spec_vals['saturated'] = \
HLP_F.get_outFillZeroSaturated(dtype)
get_outFillZeroSaturated(dtype)
self.arr.nodata = self.MetaObj.spec_vals['fill']
def calc_mean_VAA(self):
......
......@@ -41,6 +41,7 @@ from ..misc import helper_functions as HLP_F
from ..misc import path_generator as PG
from ..misc.SpatialIndexMediator import SpatialIndexMediator
from ..misc.logging import GMS_logger
from ..misc.definition_dicts import get_GMS_sensorcode, get_outFillZeroSaturated
#if socket.gethostname() == 'geoms':
......@@ -594,7 +595,7 @@ class L1B_object(L1A_object):
shift_bbl, ref_bbl = [False]*len(shift_cwl), [False]*len(ref_cwl) # bad band lists
for dic,s_r in zip([self.__dict__, ref_gmsDict], ['shift', 'ref']):
dic['GMS_identifier']['logger'] = None # set a dummy value in order to avoid Exception
sensorcode = HLP_F.get_GMS_sensorcode(dic['GMS_identifier'])
sensorcode = get_GMS_sensorcode(dic['GMS_identifier'])
if sensorcode in ['LDCM','S2A','S2B'] and '9' in dic['LayerBandsAssignment']:
locals()['%s_bbl' %s_r][dic['LayerBandsAssignment'].index('9')] = True
if sensorcode in ['S2A','S2B'] and '10' in dic['LayerBandsAssignment']:
......@@ -646,8 +647,8 @@ class L1B_object(L1A_object):
'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) for x,y in HLP_F.reorder_CornerLonLat(self.trueDataCornerLonLat)],
'nodata' : (HLP_F.get_outFillZeroSaturated(geoArr_ref .dtype)[0],
HLP_F.get_outFillZeroSaturated(geoArr_shift.dtype)[0]),
'nodata' : (get_outFillZeroSaturated(geoArr_ref .dtype)[0],
get_outFillZeroSaturated(geoArr_shift.dtype)[0]),
'ignore_errors' : True,
'v' : False,
'q' : True
......
......@@ -32,10 +32,13 @@ except ImportError:
from py_tools_ds.ptds.geo.map_info import mapinfo2geotransform
from ..config import GMS_config as CFG
from ..misc import helper_functions as HLP_F
from . import GEOPROCESSING as GEOP
from .L1B_P import L1B_object
from ..misc.definition_dicts import get_outFillZeroSaturated
from S2SCAPEM.s2_ac import AC_GMS
from S2MSI import RSImage
########################### core functions ####################################
class L1C_object(L1B_object):
......@@ -53,6 +56,7 @@ class L1C_object(L1B_object):
self.SAA_arr = None # set by self.calc_SZA_SAA_array()
self.RAA_arr = None # set by self.calc_RAA_array()
def get_lonlat_coord_array(self, subset=None):
"""Calculates pixelwise 2D-array with longitude and latitude coordinates. Supports 3 modes for subsetting:
(1) called with given subset argument if self.mask_nodata is an array: passes array subset
......@@ -64,16 +68,18 @@ class L1C_object(L1B_object):
rS, rE, cS, cE = list(GEOP.get_subsetProps_from_subsetArg(self.shape_fullArr, subset).values())[3:7]
self.log_for_fullArr_or_firstTile('Calculating lonlat array.',subset)
fillVal = HLP_F.get_outFillZeroSaturated(np.float32)[0]
fillVal = get_outFillZeroSaturated(np.float32)[0]
lonlat_arr = GEOP.get_lonlat_coord_array(self.shape_fullArr, subset[1],
mapinfo2geotransform(self.meta_odict['map info']),
self.meta_odict['coordinate system string'], self.mask_nodata, fillVal)[0]
if CFG.job.exec_mode == 'Flink' and subset[0]=='cube':
self.lonlat_arr = lonlat_arr
self.lonlat_arr = lonlat_arr
if subset[0]=='cube':
return lonlat_arr
else:
return {'desc': 'lonlat_arr', 'row_start': rS, 'row_end': rE, 'col_start': cS, 'col_end': cE,
'data': lonlat_arr}
return dict(desc='lonlat_arr',row_start=rS, row_end=rE, col_start=cS, col_end=cE, data=lonlat_arr)
def calc_acquisition_illumination_geometry(self, subset=None):
"""Calculates pixelwise arrays for viewing zenith angle (VZA), sun zenith angle (SZA),
......@@ -88,34 +94,43 @@ class L1C_object(L1B_object):
rS, rE, cS, cE = list(GEOP.get_subsetProps_from_subsetArg(self.shape_fullArr, subset).values())[3:7]
self.log_for_fullArr_or_firstTile('Calculating acquisition and illumination geometry arrays.',subset)
fillVal = HLP_F.get_outFillZeroSaturated(np.float32)[0]
fillVal = get_outFillZeroSaturated(np.float32)[0]
assert self.meta_odict, "Missing 'meta' attribute. Run self.Meta2ODict() first!"
if 'ViewingAngle_arrProv' in self.meta_odict and self.meta_odict['ViewingAngle_arrProv']:
VZA_arr = GEOP.adjust_acquisArrProv_to_shapeFullArr(self.meta_odict['ViewingAngle_arrProv'], self.shape_fullArr,
VZA_arr = GEOP.adjust_acquisArrProv_to_shapeFullArr(self.meta_odict['ViewingAngle_arrProv'],
self.shape_fullArr,
subset, bandwise=0)
else:
VZA_arr = GEOP.calc_VZA_array(self.shape_fullArr, subset[1], self.trueDataCornerPos,
float(self.meta_odict['ViewingAngle']), float(self.meta_odict['FieldOfView']),
float(self.meta_odict['ViewingAngle']),
float(self.meta_odict['FieldOfView']),
self.logger, self.mask_nodata, fillVal)
SZA_arr, SAA_arr = GEOP.calc_SZA_SAA_array(self.shape_fullArr, subset[1], self.meta_odict['AcqDate'],
self.meta_odict['AcqTime'], self.trueDataCornerPos,
self.trueDataCornerLonLat, self.meta_odict['overpass duraction sec'],
self.logger, self.mask_nodata, fillVal,
accurracy=CFG.job.SZA_SAA_calculation_accurracy,
lonlat_arr=self.lonlat_arr
if CFG.job.SZA_SAA_calculation_accurracy == 'fine' else None)
SZA_arr, SAA_arr = GEOP.calc_SZA_SAA_array(
self.shape_fullArr, subset[1],
self.meta_odict['AcqDate'],
self.meta_odict['AcqTime'],
self.trueDataCornerPos,
self.trueDataCornerLonLat,
self.meta_odict['overpass duraction sec'],
self.logger, self.mask_nodata, fillVal,
accurracy=CFG.job.SZA_SAA_calculation_accurracy,
lonlat_arr=self.lonlat_arr if CFG.job.SZA_SAA_calculation_accurracy == 'fine' else None)
if 'IncidenceAngle_arrProv' in self.meta_odict and self.meta_odict['IncidenceAngle_arrProv']:
pass # FIXME Sentinel-2s viewing azimuth array is currently not used (instead a static VAA is used)
RAA_arr = GEOP.calc_RAA_array(self.trueDataCornerLonLat, SAA_arr, self.VAA_mean, self.mask_nodata, fillVal)
if CFG.job.exec_mode == 'Flink' and subset is None:
self.VZA_arr, self.SZA_arr, self.SAA_arr, self.RAA_arr = VZA_arr, SZA_arr, SAA_arr, RAA_arr
self.VZA_arr, self.SZA_arr, self.SAA_arr, self.RAA_arr = VZA_arr, SZA_arr, SAA_arr, RAA_arr
if subset[0]=='cube':
return VZA_arr, SZA_arr, SAA_arr, RAA_arr
else:
return ({'desc': 'VZA_arr', 'row_start': rS, 'row_end': rE, 'col_start': cS, 'col_end': cE, 'data': VZA_arr},
{'desc': 'SZA_arr', 'row_start': rS, 'row_end': rE, 'col_start': cS, 'col_end': cE, 'data': SZA_arr},
{'desc': 'SAA_arr', 'row_start': rS, 'row_end': rE, 'col_start': cS, 'col_end': cE, 'data': SAA_arr},
{'desc': 'RAA_arr', 'row_start': rS, 'row_end': rE, 'col_start': cS, 'col_end': cE, 'data': RAA_arr})
return (dict(desc='VZA_arr', row_start=rS, row_end=rE, col_start=cS, col_end=cE, data=VZA_arr),
dict(desc='SZA_arr', row_start=rS, row_end=rE, col_start=cS, col_end=cE, data=SZA_arr),
dict(desc='SAA_arr', row_start=rS, row_end=rE, col_start=cS, col_end=cE, data=SAA_arr),
dict(desc='RAA_arr', row_start=rS, row_end=rE, col_start=cS, col_end=cE, data=RAA_arr))
def atm_corr(self, subset=None):
"""Performs an atmospheric correction and returns atmospherically corrected reflectance data."""
......@@ -129,6 +144,7 @@ class L1C_object(L1B_object):
self.delete_acquisition_illumination_geometry()
self.lonlat_arr = None
def delete_acquisition_illumination_geometry(self):
self.VZA_arr = None # not needed anymore
self.SZA_arr = None # not needed anymore
......@@ -155,15 +171,17 @@ class AtmCorr(object):
self._metadata = {}
self._nodata = {}
self._yesdata = {}
self._band_spatial_sampling = {}
# assertions
scene_IDs = [obj.scene_ID for obj in L1C_objs]
assert len(list(set(scene_IDs)))==1, \
"Input GMS objects for 'AtmCorr' must all belong to the same scene ID!. Received %s." %scene_IDs
self.inObjs = L1C_objs
self.results = None # direct output of external atmCorr module (set by run_atmospheric_correction)
self.outObjs = None # atmospherically corrected L1C objects
self.inObjs = L1C_objs
self.ac_input = {} # set by self.run_atmospheric_correction()
self.results = None # direct output of external atmCorr module (set by run_atmospheric_correction)
self.outObjs = None # atmospherically corrected L1C objects
@property
......@@ -204,13 +222,15 @@ class AtmCorr(object):
if not self._data:
for inObj in self.inObjs:
for bandIdx, bandN in enumerate(inObj.LayerBandsAssignment):
for bandN, bandIdx in inObj.arr.bandnames.items():
if bandN not in self._data:
# TODO adjust band name to match B0X format
self._data[bandN] = inObj.arr[:,:,bandIdx] # TODO arr should also support band names from LayerBandsAssignment
arr2pass = inObj.arr[:,:,bandIdx].astype(np.float32) # conversion to np.float16 will convert -9999 to -10000
arr2pass[arr2pass==inObj.arr.nodata] = np.nan # set nodata values to np.nan
self._data[bandN] = (arr2pass/inObj.meta_odict['ScaleFactor']).astype(np.float16)
else:
inObj.logger.warning(
"Band '%s' cannot be included into atmospheric correction because it exists multiple times.")
inObj.logger.warning("Band '%s' cannot be included into atmospheric correction because it "
"exists multiple times.")
return self._data
......@@ -222,6 +242,63 @@ class AtmCorr(object):
self._data = data_dict
@property
def nodata(self):
"""
:return:
___ attribute: nodata, type:<class 'dict'>
______ key:60.0, value_type:<class 'numpy.ndarray'>, repr: [[ TrueTrueTrue ..., [..] False False False]]
______ key:10.0, value_type:<class 'numpy.ndarray'>, repr: [[ TrueTrueTrue ..., [..] False False False]]
______ key:20.0, value_type:<class 'numpy.ndarray'>, repr: [[ TrueTrueTrue ..., [..] False False False]]
"""
if not self._nodata:
for inObj in self.inObjs:
self._nodata[inObj.arr.xgsd] = ~inObj.arr.mask_nodata[:]
return self._nodata
@property
def tile_name(self):
"""Returns S2A tile name
:return: e.g.
'32UMA'
"""
return '' # FIXME
@property
def band_spatial_sampling(self):
"""
:return: e.g.
{'B01': 60.0,
'B02': 10.0,
'B03': 10.0,
'B04': 10.0,
'B05': 20.0,
'B06': 20.0,
'B07': 20.0,
'B08': 10.0,
'B09': 60.0,
'B10': 60.0,
'B11': 20.0,
'B12': 20.0,
'B8A': 20.0}
"""
if not self._band_spatial_sampling:
for inObj in self.inObjs:
for bandN in inObj.arr.bandnames:
if bandN not in self._band_spatial_sampling:
self._band_spatial_sampling[bandN] = inObj.arr.xgsd
return self._band_spatial_sampling
@property
def metadata(self):
"""
......@@ -254,75 +331,190 @@ class AtmCorr(object):
"""
if not self._metadata:
# set corner coordinates and dims
self._metadata['spatial_samplings'] = {}
for inObj in self.inObjs:
# get GSD
if inObj.arr.xgsd != inObj.arr.ygsd:
warnings.warn("X/Y GSD is not equal for entity ID %s" % inObj.entity_ID +
(' (%s)' % inObj.subsystem if inObj.subsystem else '') +
'Using X-GSD as key for spatial sampling dictionary.')
gsd = inObj.arr.xgsd
# set spatial information
self._metadata['spatial_samplings'][gsd] = dict(
ULX = inObj.arr.box.boxMapYX[0][1],
ULY = inObj.arr.box.boxMapYX[0][0],
XDIM = inObj.arr.xgsd,
YDIM = inObj.arr.ygsd,
NROWS = inObj.arr.rows,
NCOLS = inObj.arr.cols)
# set sensing time
self._metadata['SENSING_TIME'] = self.inObjs[0].acquisition_date
self._metadata['U'] = self.inObjs[0].meta_odict['EarthSunDist']
self._metadata['SENSING_TIME'] = self.inObjs[0].acquisition_date
self._metadata['viewing_zenith'] = self._meta_get_viewing_zenith()
self._metadata['viewing_azimuth'] = self._meta_get_viewing_azimuth()
self._metadata['relative_viewing_azimuth'] = self._meta_get_relative_viewing_azimuth()
self._metadata['sun_mean_azimuth'] = self.inObjs[0].meta_odict['SunAzimuth']
self._metadata['sun_mean_zenith'] = 90-self.inObjs[0].meta_odict['SunElevation']
self._metadata['solar_irradiance'] = self._meta_get_solar_irradiance()
self._metadata['aux_data'] = {}
self._metadata['spatial_samplings'] = self._meta_get_spatial_samplings()
return self._metadata
@property
def nodata(self):
def _meta_get_spatial_samplings(self):
"""
:return:
___ attribute: nodata, type:<class 'dict'>
______ key:60.0, value_type:<class 'numpy.ndarray'>, repr: [[ TrueTrueTrue ..., [..] False False False]]
______ key:10.0, value_type:<class 'numpy.ndarray'>, repr: [[ TrueTrueTrue ..., [..] False False False]]
______ key:20.0, value_type:<class 'numpy.ndarray'>, repr: [[ TrueTrueTrue ..., [..] False False False]]
{10.0: {'NCOLS': 10980,
'NROWS': 10980,
'ULX': 499980.0,
'ULY': 5800020.0,
'XDIM': 10.0,
'YDIM': -10.0},
20.0: {'NCOLS': 5490,
'NROWS': 5490,
'ULX': 499980.0,
'ULY': 5800020.0,
'XDIM': 20.0,
'YDIM': -20.0},
60.0: {'NCOLS': 1830,
'NROWS': 1830,
'ULX': 499980.0,
'ULY': 5800020.0,
'XDIM': 60.0,
'YDIM': -60.0}}
"""
# set corner coordinates and dims
spatial_samplings = {}
for inObj in self.inObjs:
self._nodata[inObj.arr.xgsd] = ~inObj.arr.mask_nodata[:]
return self._nodata
# validate GSD
if inObj.arr.xgsd != inObj.arr.ygsd:
warnings.warn("X/Y GSD is not equal for entity ID %s" % inObj.entity_ID +
(' (%s)' % inObj.subsystem if inObj.subsystem else '') +
'Using X-GSD as key for spatial sampling dictionary.')
# set spatial information
spatial_samplings[inObj.arr.xgsd] = dict(
ULX=inObj.arr.box.boxMapYX[0][1],
ULY=inObj.arr.box.boxMapYX[0][0],
XDIM=inObj.arr.xgsd,
YDIM=-inObj.arr.ygsd,
NROWS=inObj.arr.rows,
NCOLS=inObj.arr.cols)
@property
def yesdata(self):
return spatial_samplings
def _meta_get_solar_irradiance(self):
"""
:return:
___ atribute: yesdata, type:<class 'dict'>
______ key:20.0, value_type:<class 'numpy.ndarray'>, repr: [[False False False [..] e ...,TrueTrueTrue]]
______ key:10.0, value_type:<class 'numpy.ndarray'>, repr: [[False False False [..] e ...,TrueTrueTrue]]
______ key:60.0, value_type:<class 'numpy.ndarray'>, repr: [[False False False [..] e ...,TrueTrueTrue]]
{'B01': 1913.57,
'B02': 1941.63,
'B03': 1822.61,
'B04': 1512.79,
'B05': 1425.56,
'B06': 1288.32,
'B07': 1163.19,
'B08': 1036.39,
'B09': 813.04,
'B10': 367.15,
'B11': 245.59,
'B12': 85.25,
'B8A': 955.19}
"""
solar_irradiance = {}
for inObj in self.inObjs:
for bandN, bandIdx in inObj.arr.bandnames.items():
if bandN not in solar_irradiance:
solar_irradiance[bandN] = inObj.meta_odict['SolIrradiance'][bandIdx]
return solar_irradiance
def _meta_get_viewing_zenith(self):
"""
:return: {B10:ndarray(dtype=float16),[...],B09:ndarray(dtype=float16)}
"""
viewing_zenith = {}
for inObj in self.inObjs:
if inObj.VZA_arr is None:
inObj.calc_acquisition_illumination_geometry()
for bandN, bandIdx in inObj.arr.bandnames.items():
if bandN not in viewing_zenith:
viewing_zenith[bandN] = inObj.VZA_arr[:, :, bandIdx] if inObj.VZA_arr.ndim == 3 else inObj.VZA_arr
return viewing_zenith
def _meta_get_viewing_azimuth(self):
"""
:return: {B10:ndarray(dtype=float16),[...],B09:ndarray(dtype=float16)}
"""