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

first (not yet working) version of wrapper for atmospheric correction

algorithms.__init__:
- added some imports
algorithms.gms_object:
 - removed deprecated code
algorithms.L1C_P:
- removed deprecated code
- added class 'AtmCorr': a wrapper class for atmospheric correction
misc.database_tools:
- fixed some broken type hints
misc.helper_functions
- fixed broken type hint
- replaced deprecated function call
processing.pipeline:
- fixed some broken type hints
- updated __version__
parent 428ee490
......@@ -15,7 +15,7 @@ from . import config
from .processing.process_controller import process_controller
__version__ = '20170105.01'
__version__ = '20170106.01'
__author__ = 'Daniel Scheffler'
__all__ = ['algorithms',
'io',
......
......@@ -21,6 +21,8 @@
###############################################################################
########################### Library import ####################################
import warnings
import numpy as np
try:
from osgeo import osr
......@@ -33,7 +35,7 @@ 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 S2SCAPEM.s2_ac import AC_GMS
########################### core functions ####################################
class L1C_object(L1B_object):
......@@ -62,12 +64,6 @@ 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)
#if hasattr(self, 'mask_nodata') and isinstance(self.mask_nodata, np.ndarray):
# mask_1bit_temp = self.mask_nodata[rS:rE + 1, cS:cE + 1]
#else:
# path_masks = PG.path_generator(self.__dict__.copy(),proc_level='L1B').get_path_maskdata()
# mask_1bit_temp = INP_R.read_mask_subset(path_masks, 'mask_nodata', self.logger, subset)
fillVal = HLP_F.get_outFillZeroSaturated(np.float32)[0]
lonlat_arr = GEOP.get_lonlat_coord_array(self.shape_fullArr, subset[1],
mapinfo2geotransform(self.meta_odict['map info']),
......@@ -92,12 +88,6 @@ 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)
#if hasattr(self, 'mask_nodata') and isinstance(self.mask_nodata, np.ndarray):
# mask_1bit_temp = self.mask_nodata[rS:rE + 1, cS:cE + 1]
#else:
# path_masks = PG.path_generator(self.__dict__.copy(),proc_level='L1B').get_path_maskdata()
# mask_1bit_temp = INP_R.read_mask_subset(path_masks, 'mask_nodata', self.logger, subset)
fillVal = HLP_F.get_outFillZeroSaturated(np.float32)[0]
assert self.meta_odict, "Missing 'meta' attribute. Run self.Meta2ODict() first!"
......@@ -145,3 +135,194 @@ class L1C_object(L1B_object):
self.SAA_arr = None # not needed anymore
self.RAA_arr = None # not needed anymore
class AtmCorr(object):
def __init__(self, *L1C_objs):
"""Wrapper for atmospheric correction, written by Andre Hollstein, GFZ Potsdam.
Creates the input arguments for atmospheric correction from one or multiple L1C_object instance(s) belonging to
the same scene ID, performs the atmospheric correction and returns the atmospherically corrected L1C object(s).
:param L1C_objs: one or more instances of L1C_object belonging to the same scene ID
"""
L1C_objs = L1C_objs if isinstance(L1C_objs, tuple) else (L1C_objs,)
# hidden attributes
self._GSDs = []
self._data = {}
self._metadata = {}
self._nodata = {}
self._yesdata = {}
# 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
@property
def GSDs(self):
"""
Returns a list of spatial samplings within the input GMS objects, e.g. [10,20,60].
"""
for obj in self.inObjs:
if obj.arr.xgsd != obj.arr.ygsd:
warnings.warn("X/Y GSD is not equal for entity ID %s" % obj.entity_ID +
(' (%s)'%obj.subsystem if obj.subsystem else '') +
'Using X-GSD as key for spatial sampling dictionary.')
self._GSDs.append(obj.arr.xgsd)
return self._GSDs
@property
def data(self):
"""
:return:
___ attribute: data, type:<class 'dict'>
______ key:B05, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 085998540.0803833 ]]
______ key:B01, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 131225590.13208008]]
______ key:B06, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] .14965820.13977051]]
______ key:B11, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] .11492920.10192871]]
______ key:B02, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 104187010.10308838]]
______ key:B10, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 013099670.01300049]]
______ key:B08, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] .16857910.15783691]]
______ key:B04, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 065490720.06228638]]
______ key:B03, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 082702640.08148193]]
______ key:B12, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 068420410.06060791]]
______ key:B8A, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 192138670.17553711]]
______ key:B09, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] .09600830.09887695]]
______ key:B07, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 173339840.15600586]]
"""
if not self._data:
for inObj in self.inObjs:
for bandIdx, bandN in enumerate(inObj.LayerBandsAssignment):
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
else:
inObj.logger.warning(
"Band '%s' cannot be included into atmospheric correction because it exists multiple times.")
return self._data
@data.setter
def data(self, data_dict):
assert isinstance(data_dict, dict), \
"'data' can only be set to a dictionary with band names as keys and numpy arrays as values."
self._data = data_dict
@property
def metadata(self):
"""
:return:
___ attribute: metadata, type:<class 'dict'>
______ key:spatial_samplings
_________ key:60.0
____________ key:ULY, value_type:<class 'int'>, repr: 4900020
____________ key:NCOLS, value_type:<class 'int'>, repr: 1830
____________ key:XDIM, value_type:<class 'int'>, repr: 60
____________ key:ULX, value_type:<class 'int'>, repr: 600000
____________ key:NROWS, value_type:<class 'int'>, repr: 1830
____________ key:YDIM, value_type:<class 'int'>, repr: -60
_________ key:10.0
____________ key:ULY, value_type:<class 'int'>, repr: 4900020
____________ key:NCOLS, value_type:<class 'int'>, repr: 10980
____________ key:XDIM, value_type:<class 'int'>, repr: 10
____________ key:ULX, value_type:<class 'int'>, repr: 600000
____________ key:NROWS, value_type:<class 'int'>, repr: 10980
____________ key:YDIM, value_type:<class 'int'>, repr: -10
_________ key:20.0
____________ key:ULY, value_type:<class 'int'>, repr: 4900020
____________ key:NCOLS, value_type:<class 'int'>, repr: 5490
____________ key:XDIM, value_type:<class 'int'>, repr: 20
____________ key:ULX, value_type:<class 'int'>, repr: 600000
____________ key:NROWS, value_type:<class 'int'>, repr: 5490
____________ key:YDIM, value_type:<class 'int'>, repr: -20
______ key:SENSING_TIME, value_type:<class 'datetime.datetime'>, repr: 2016-03-26 10:34:06.538000+00:00
"""
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
return self._metadata
@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]]
"""
for inObj in self.inObjs:
self._nodata[inObj.arr.xgsd] = ~inObj.arr.mask_nodata[:]
return self._nodata
@property
def yesdata(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]]
"""
for inObj in self.inObjs:
self._yesdata[inObj.arr.xgsd] = inObj.arr.mask_nodata[:]
return self._yesdata
def run_atmospheric_correction(self):
inArgs = self.data, self.metadata, self.nodata, self.yesdata
self.results = AC_GMS(inArgs)
return
def join_results_to_inObjs(self):
self.outObjs = self.results # TODO
return self.outObjs
......@@ -1440,6 +1440,7 @@ class METADATA(object):
def get_EarthSunDistance(self,acqDate):
"""Get earth sun distance (requires file of pre calculated earth sun distance per day)
:param acqDate:
"""
if not os.path.exists(CFG.job.path_earthSunDist):
......@@ -1459,7 +1460,9 @@ class METADATA(object):
return float(EA_dist_dict[acqDate])
def calc_center_acquisition_time(self,trueDataCornerLonLat,logger):
"""Calculates a mission center acquision time using acquisition date, true date corner coordinates and solar azimuth.
"""
Calculates a mission center acquision time using acquisition date, true date corner coordinates and solar azimuth.
:param trueDataCornerLonLat:
:param logger:
"""
......@@ -1482,6 +1485,7 @@ class METADATA(object):
def get_overpassDuration_SceneLength(self, trueDataCornerLonLat, trueDataCornerPos, shape_fullArr,logger):
"""Calculates duration of image acquisition in seconds.
:param trueDataCornerLonLat:
:param trueDataCornerPos:
:param shape_fullArr:
......@@ -1516,6 +1520,7 @@ class METADATA(object):
def get_LayerBandsAssignment(GMS_identifier, nBands=None, ignore_usecase=False):
"""Returns LayerBandsAssignment corresponding to given satellite, sensor and subsystem and with respect to
CFG.usecase.sort_bands_by_cwl, CFG.usecase.skip_thermal and CFG.usecase.skip_pan.
:param GMS_identifier: <dict>, derived from self.get_GMS_identifier()
:param nBands: should be specified if number of bands differs from standard
(e.g. SPOT data containing no PAN)
......
......@@ -7,29 +7,24 @@
#
###############################################################################
# FIXME comment that in as soon as config is no longer part of builtins
# from . import GEOPROCESSING
# from . import METADATA
# from . import gms_cloud_classifier
# from . import L0A_P
# from . import L0B_P
# from . import L1A_P
# from . import L1B_P
# from . import L1C_P
# from . import L2A_P
# from . import L2B_P
# from . import L2C_P
#
# __all__=['GEOPROCESSING',
# 'METADATA',
# 'gms_cloud_classifier',
# 'L0A_P',
# 'L0B_P',
# 'L1A_P',
# 'L1B_P',
# 'L1C_P',
# 'L2A_P',
# 'L2B_P',
# 'L2C_P']
from . import GEOPROCESSING
from . import METADATA
from . import gms_cloud_classifier
from . import L1A_P
from . import L1B_P
from . import L1C_P
from . import L2A_P
from . import L2B_P
from . import L2C_P
__all__=['GEOPROCESSING',
'METADATA',
'gms_cloud_classifier',
'L1A_P',
'L1B_P',
'L1C_P',
'L2A_P',
'L2B_P',
'L2C_P']
__author__='Daniel Scheffler'
......@@ -66,7 +66,7 @@ class GMS_object(object):
self.sensor = ''
self.subsystem = ''
self.sensormode = ''
self.acquisition_date = None
self.acquisition_date = None # also includes time, set by from_disk() and within L1A_P
self.entity_ID = ''
self.scene_ID = -9999
self.filename = ''
......@@ -494,18 +494,6 @@ class GMS_object(object):
self.arr = self.pathGen.get_path_imagedata()
self.masks = self.pathGen.get_path_maskdata() # self.mask_nodata and self.mask_clouds are auto-synchronized via self.masks (see their getters)
#if CFG.job.exec_mode=='Flink':
# # TODO use GeoArray read functions (otherwise map info is not used)
# self.arr = INP_R.read_ENVIfile(path_arr,self.arr_shape,self.arr_pos,self.logger,q=1)
# self.mask_nodata = INP_R.read_mask_subset(path_masks, 'mask_nodata', self.logger, tuple_GMS_subset[1])
# self.mask_clouds = INP_R.read_mask_subset(path_masks, 'mask_clouds', self.logger, tuple_GMS_subset[1])
# self.log_for_fullArr_or_firstTile('Reading file %s as tiles...' %self.baseN \
# if self.arr_pos else 'Reading file %s...' %self.baseN)
# #self.masks is only needed by Output writer to masks combined -> generated there and on demand
#else: # CFG.job.exec_mode=='Python'
# self.arr = path_arr
# self.masks = path_masks # self.mask_nodata and self.mask_clouds are auto-synchronized via self.masks (see their getters)
return copy.copy(self)
......
......@@ -94,7 +94,7 @@ def get_info_from_SQLdb(path_db,tablename,vals2return,cond_dict,records2fetch=0)
def get_postgreSQL_value(value):
# type: (str,Any) -> str
# type: (Any) -> str
"""Converts Python variable to a postgreSQL value respecting postgreSQL type casts.
The resulting value can be directly inserted into a postgreSQL query."""
......@@ -169,7 +169,7 @@ def get_info_from_postgreSQLdb(conn_params,tablename,vals2return,cond_dict=None,
def update_records_in_postgreSQLdb(conn_params, tablename, vals2update_dict, cond_dict=None, timeout=15000):
# type: (str, str, dict, dict, int)
# type: (str, str, dict, dict, int) -> None
"""Queries a postgreSQL database for the given parameters and updates the given columns of the query result.
:param conn_params: <str> connection parameters as provided by CFG.job.conn_params
......@@ -200,7 +200,7 @@ def update_records_in_postgreSQLdb(conn_params, tablename, vals2update_dict, con
def append_item_to_arrayCol_in_postgreSQLdb(conn_params, tablename, vals2append_dict, cond_dict=None, timeout=15000):
# type: (str, str, dict, dict, int)
# type: (str, str, dict, dict, int) -> None
"""Queries a postgreSQL database for the given parameters
and appends the given value to the specified column of the query result.
......@@ -238,7 +238,7 @@ def append_item_to_arrayCol_in_postgreSQLdb(conn_params, tablename, vals2append_
def remove_item_from_arrayCol_in_postgreSQLdb(conn_params, tablename, vals2remove_dict, cond_dict=None, timeout=15000):
# type: (str, str, dict, dict, int)
# type: (str, str, dict, dict, int) -> None
"""Queries a postgreSQL database for the given parameters
and removes the given value from the specified column of the query result.
......@@ -276,7 +276,7 @@ def remove_item_from_arrayCol_in_postgreSQLdb(conn_params, tablename, vals2remov
def create_record_in_postgreSQLdb(conn_params, tablename, vals2write_dict, timeout=15000):
# type: (str, str, dict, dict, int)
# type: (str, str, dict, dict, int) -> int
"""Creates a single new record in a postgreSQL database and pupulates its columns with the given values.
:param conn_params: <str> connection parameters as provided by CFG.job.conn_params
......
......@@ -57,7 +57,7 @@ from ..algorithms.L2A_P import L2A_object
from ..algorithms.L2B_P import L2B_object
from ..algorithms.L2C_P import L2C_object
from py_tools_ds.ptds.geo.coord_trafo import mapXY2imXY, reproject_shapelyPoly
from py_tools_ds.ptds.geo.coord_trafo import mapXY2imXY, reproject_shapelyGeometry
from py_tools_ds.ptds.geo.coord_calc import corner_coord_to_minmax
......@@ -406,7 +406,7 @@ def get_image_tileborders(target_tileShape, target_tileSize, path_GMS_file=None,
def cut_GMS_obj_into_blocks(tuple__In_obj__blocksize_RowsCols): # FIXME should better call get_subset_GMS_obj (also takes care to tile map infos)
# type: (tuple) -> L1A_object
# type: (tuple) -> list
"""Cut a GMS object into tiles with respect to raster attributes as well as scene wide attributes.
:param tuple__In_obj__blocksize_RowsCols: a tuple with GMS_obj as first and [rows,cols] as second element"""
......@@ -680,7 +680,7 @@ def get_arrSubsetBounds_from_shapelyPolyLonLat(arr_shape, shpPolyLonLat, im_gt,
:param ensure_valid_coords: <bool> whether to ensure that the returned values are all inside the original
image bounding box
"""
shpPolyImPrj = reproject_shapelyPoly(shpPolyLonLat,im_prj)
shpPolyImPrj = reproject_shapelyGeometry(shpPolyLonLat, 4326, im_prj)
imCoordsXY = get_imageCoords_from_shapelyPoly(shpPolyImPrj,im_gt)
bounds = corner_coord_to_minmax(imCoordsXY)
outbounds = get_valid_arrSubsetBounds(arr_shape, bounds, buffer=pixbuffer) if ensure_valid_coords else bounds
......
......@@ -50,7 +50,7 @@ def L1A_map(dataset_dict): #map (scene-wise parallelization)
@HLP_F.log_uncaught_exceptions
def L1A_map_1(dataset_dict, block_size=None): #map (scene-wise parallelization)
# type: (dict) -> L1A_P.L1A_object
# type: (dict) -> list
L1A_obj = L1A_P.L1A_object(image_type = 'RSD',
satellite = dataset_dict['satellite'],
......@@ -131,7 +131,7 @@ def L1C_map(L1B_obj):
@HLP_F.log_uncaught_exceptions
def L2A_map(L1C_obj, block_size=None):
# type: (L1C_P.L1C_object, tuple) -> L2A_P.L2A_object
# type: (L1C_P.L1C_object, tuple) -> list
L2A_obj = L2A_P.L2A_object(L1C_obj)
L2A_obj.correct_spatial_shifts()
L2A_obj.calc_mask_nodata(overwrite=True) # update no data mask
......
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