Commit 65de3567 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Revised L1B_P.Scene_finder() and L1B_P.L1B_object.get_spatial_reference_scene()

parent d81b937f
Pipeline #1252 failed with stage
in 8 minutes and 51 seconds
...@@ -39,13 +39,14 @@ __author__ = 'Daniel Scheffler' ...@@ -39,13 +39,14 @@ __author__ = 'Daniel Scheffler'
class Scene_finder(object): class Scene_finder(object):
def __init__(self, src_boundsLonLat, src_AcqDate, src_prj, src_footprint_poly, min_overlap=20, def __init__(self, src_boundsLonLat, src_AcqDate, src_prj, src_footprint_poly, sceneID_excluded=None,
min_cloudcov=0, max_cloudcov=20, plusminus_days=30, plusminus_years=10): min_overlap=20, min_cloudcov=0, max_cloudcov=20, plusminus_days=30, plusminus_years=10):
self.boundsLonLat = src_boundsLonLat self.boundsLonLat = src_boundsLonLat
self.src_AcqDate = src_AcqDate self.src_AcqDate = src_AcqDate
self.src_prj = src_prj self.src_prj = src_prj
self.src_footprint_poly = src_footprint_poly self.src_footprint_poly = src_footprint_poly
self.sceneID_excluded = sceneID_excluded
self.min_overlap = min_overlap self.min_overlap = min_overlap
self.min_cloudcov = min_cloudcov self.min_cloudcov = min_cloudcov
self.max_cloudcov = max_cloudcov self.max_cloudcov = max_cloudcov
...@@ -99,106 +100,103 @@ class Scene_finder(object): ...@@ -99,106 +100,103 @@ class Scene_finder(object):
GDF['polyUTM'] = list(GDF['polyLonLat'].map(LonLat2UTM)) GDF['polyUTM'] = list(GDF['polyLonLat'].map(LonLat2UTM))
self.GDF_ref_scenes = GDF self.GDF_ref_scenes = GDF
def filter_possib_ref_scenes(self): def _collect_refscene_metadata(self):
"""Filter possible scenes by running all filter functions.""" """Collect some reference scene metadata needed for later filtering."""
self._filter_by_overlap() GDF = self.GDF_ref_scenes
self._filter_by_proc_status()
self._filter_by_dataset_existance()
self._filter_by_entity_ID_availability()
self._filter_by_projection()
def choose_ref_scene(self):
"""Choose reference scene with minimum cloud cover and maximum overlap."""
if not self.GDF_ref_scenes.empty:
GDF = self.GDF_ref_scenes
GDF = GDF[GDF['cloudcover'] == GDF['cloudcover'].min()]
GDF = GDF[GDF['overlap percentage'] == GDF['overlap percentage'].max()]
if not GDF.empty: # get overlap parameter
GDF_res = GDF.iloc[0] def get_OL_prms(poly): return get_overlap_polygon(poly, self.src_footprint_poly)
return ref_Scene(GDF_res)
else: GDF['overlapParams'] = list(GDF['polyLonLat'].map(get_OL_prms))
return None GDF['overlap area'] = list(GDF['overlapParams'].map(lambda OL_prms: OL_prms['overlap area']))
GDF['overlap percentage'] = list(GDF['overlapParams'].map(lambda OL_prms: OL_prms['overlap percentage']))
GDF['overlap poly'] = list(GDF['overlapParams'].map(lambda OL_prms: OL_prms['overlap poly']))
del GDF['overlapParams']
# get processing level of reference scenes
procL = GeoDataFrame(
DB_T.get_info_from_postgreSQLdb(CFG.job.conn_database, 'scenes_proc', ['sceneid', 'proc_level'],
{'sceneid': list(GDF.sceneid)}), columns=['sceneid', 'proc_level'])
GDF = GDF.merge(procL, on='sceneid', how='left')
GDF = GDF.where(GDF.notnull(), None) # replace NaN values with None
# get path of binary file
def get_path_binary(GDF_row):
return PG.path_generator(scene_ID=GDF_row['sceneid'], proc_level=GDF_row['proc_level']) \
.get_path_imagedata() if GDF_row['proc_level'] else None
GDF['path_ref'] = GDF.apply(lambda GDF_row: get_path_binary(GDF_row), axis=1)
GDF['refDs_exists'] = list(GDF['path_ref'].map(lambda p: os.path.exists(p) if p else False))
# check if a proper entity ID can be gathered from database
eID = GeoDataFrame(DB_T.get_info_from_postgreSQLdb(CFG.job.conn_database, 'scenes', ['id', 'entityid'],
{'id': list(GDF.sceneid)}), columns=['sceneid', 'entityid'])
GDF = GDF.merge(eID, on='sceneid', how='left')
self.GDF_ref_scenes = GDF.where(GDF.notnull(), None)
def _filter_excluded_sceneID(self):
"""Filter reference scene with the same scene ID like the target scene."""
GDF = self.GDF_ref_scenes
if not GDF.empty:
self.GDF_ref_scenes = GDF.loc[GDF['sceneid'] != self.sceneID_excluded]
def _filter_by_overlap(self): def _filter_by_overlap(self):
"""Filter all scenes with less spatial overlap than self.min_overlap.""" """Filter all scenes with less spatial overlap than self.min_overlap."""
GDF = self.GDF_ref_scenes GDF = self.GDF_ref_scenes
if not GDF.empty: if not GDF.empty:
# get overlap parameter
def get_OL_prms(poly): return get_overlap_polygon(poly, self.src_footprint_poly)
GDF['overlapParams'] = list(GDF['polyLonLat'].map(get_OL_prms))
GDF['overlap area'] = list(GDF['overlapParams'].map(lambda OL_prms: OL_prms['overlap area']))
GDF['overlap percentage'] = list(GDF['overlapParams'].map(lambda OL_prms: OL_prms['overlap percentage']))
GDF['overlap poly'] = list(GDF['overlapParams'].map(lambda OL_prms: OL_prms['overlap poly']))
del GDF['overlapParams']
# filter scenes out that have too less overlap
self.GDF_ref_scenes = GDF.loc[GDF['overlap percentage'] >= self.min_overlap] self.GDF_ref_scenes = GDF.loc[GDF['overlap percentage'] >= self.min_overlap]
def _filter_by_proc_status(self): def _filter_by_proc_status(self):
"""Filter all scenes that have not been processed before according to proc. status (at least L1A is needed).""" """Filter all scenes that have not been processed before according to proc. status (at least L1A is needed)."""
GDF = self.GDF_ref_scenes GDF = self.GDF_ref_scenes
if not GDF.empty: if not GDF.empty:
# get processing level of refernce scenes
def query_procL(sceneID):
return DB_T.get_info_from_postgreSQLdb(CFG.job.conn_database, 'scenes_proc', ['proc_level'],
{'sceneid': sceneID})
GDF['temp_queryRes'] = list(GDF['sceneid'].map(query_procL))
GDF['proc_level'] = list(GDF['temp_queryRes'].map(lambda queryRes: queryRes[0][0] if queryRes else None))
GDF.drop('temp_queryRes', axis=1, inplace=True)
# filter scenes out that have not yet been processed
self.GDF_ref_scenes = GDF[GDF['proc_level'].notnull()] self.GDF_ref_scenes = GDF[GDF['proc_level'].notnull()]
def _filter_by_dataset_existance(self): def _filter_by_dataset_existance(self):
"""Filter all scenes where no processed data can be found on disk.""" """Filter all scenes where no processed data can be found on fileserver."""
GDF = self.GDF_ref_scenes GDF = self.GDF_ref_scenes
if not GDF.empty: if not GDF.empty:
# get path of binary file and check if the corresponding dataset exists
GDF = self.GDF_ref_scenes
def get_path_binary(GDF_row):
return PG.path_generator(scene_ID=GDF_row['sceneid'], proc_level=GDF_row['proc_level'])\
.get_path_imagedata()
def check_exists(path): return os.path.exists(path)
GDF['path_ref'] = GDF.apply(lambda GDF_row: get_path_binary(GDF_row), axis=1)
GDF['refDs_exists'] = list(GDF['path_ref'].map(check_exists))
# filter scenes out where the corresponding dataset does not exist on fileserver
self.GDF_ref_scenes = GDF[GDF['refDs_exists']] self.GDF_ref_scenes = GDF[GDF['refDs_exists']]
def _filter_by_entity_ID_availability(self): def _filter_by_entity_ID_availability(self):
"""Filter all scenes where no proper entity ID can be found in the database.""" """Filter all scenes where no proper entity ID can be found in the database (database errors)."""
GDF = self.GDF_ref_scenes GDF = self.GDF_ref_scenes
if not GDF.empty: if not GDF.empty:
# check if a proper entity ID can be gathered from database
def query_eID(sceneID):
return DB_T.get_info_from_postgreSQLdb(CFG.job.conn_database, 'scenes', ['entityid'],
{'id': sceneID}, records2fetch=1)
GDF['temp_queryRes'] = list(GDF['sceneid'].map(query_eID))
GDF['entityid'] = list(GDF['temp_queryRes'].map(lambda queryRes: queryRes[0][0] if queryRes else None))
GDF.drop('temp_queryRes', axis=1, inplace=True)
# filter scenes out that have no entity ID (database errors)
self.GDF_ref_scenes = GDF[GDF['entityid'].notnull()] self.GDF_ref_scenes = GDF[GDF['entityid'].notnull()]
def _filter_by_projection(self): def _filter_by_projection(self):
"""Filter all scenes that have a different projection than the target image.""" """Filter all scenes that have a different projection than the target image."""
GDF = self.GDF_ref_scenes GDF = self.GDF_ref_scenes[self.GDF_ref_scenes.refDs_exists]
if not GDF.empty: if not GDF.empty:
# compare projections of target and reference image # compare projections of target and reference image
from ..io.input_reader import read_ENVIhdr_to_dict GDF['prj_equal'] = \
list(GDF['path_ref'].map(lambda path_ref: prj_equal(self.src_prj, GeoArray(path_ref).prj)))
self.GDF_ref_scenes = GDF[GDF['prj_equal']]
def get_prj(path_binary): def choose_ref_scene(self):
return read_ENVIhdr_to_dict(os.path.splitext(path_binary)[0] + '.hdr')['coordinate system string'] """Choose reference scene with minimum cloud cover and maximum overlap."""
# First, collect some relavant reference scene metadata
self._collect_refscene_metadata()
def is_prj_equal(path_binary): return prj_equal(self.src_prj, get_prj(path_binary)) # Filter possible scenes by running all filter functions
GDF['prj_equal'] = list(GDF['path_ref'].map(is_prj_equal)) self._filter_excluded_sceneID()
self._filter_by_overlap()
self._filter_by_proc_status()
self._filter_by_dataset_existance()
self._filter_by_entity_ID_availability()
self._filter_by_projection()
# filter scenes out that have a different projection # Choose the reference scene out of the filtered DataFrame
self.GDF_ref_scenes = GDF[GDF['prj_equal']] if not self.GDF_ref_scenes.empty:
GDF = self.GDF_ref_scenes
GDF = GDF[GDF['cloudcover'] == GDF['cloudcover'].min()]
GDF = GDF[GDF['overlap percentage'] == GDF['overlap percentage'].max()]
if not GDF.empty:
GDF_res = GDF.iloc[0]
return ref_Scene(GDF_res)
else:
return None
class ref_Scene: class ref_Scene:
def __init__(self, GDF_record): def __init__(self, GDF_record):
...@@ -245,7 +243,7 @@ class L1B_object(L1A_object): ...@@ -245,7 +243,7 @@ class L1B_object(L1A_object):
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.meta_odict['coordinate system string'],
footprint_poly, 20, 0, 20, 30, 10) footprint_poly, self.scene_ID, 20, 0, 20, 30, 10)
# run spatial query # run spatial query
self.logger.info('Querying database in order to find a suitable reference scene for co-registration.') self.logger.info('Querying database in order to find a suitable reference scene for co-registration.')
...@@ -255,21 +253,17 @@ class L1B_object(L1A_object): ...@@ -255,21 +253,17 @@ class L1B_object(L1A_object):
else: else:
self.logger.warning('Spatial query returned no matches. Coregistration impossible.') self.logger.warning('Spatial query returned no matches. Coregistration impossible.')
self.spatRef_available = False self.spatRef_available = False
return None
# filter results # try to get a spatial reference scene by applying some filter criteria
RSF.filter_possib_ref_scenes() self.spatRef_scene = RSF.choose_ref_scene()
if RSF.GDF_ref_scenes.empty: if self.spatRef_scene:
self.spatRef_available = True
self.logger.info('Found a suitable reference image for coregistration: scene ID %s (entity ID %s).'
% (self.spatRef_scene.scene_ID, self.spatRef_scene.entity_ID))
else:
self.spatRef_available = False
self.logger.warning('No scene fulfills all requirements to serve as spatial reference for scene %s ' self.logger.warning('No scene fulfills all requirements to serve as spatial reference for scene %s '
'(entity ID %s). Coregistration impossible.' % (self.scene_ID, self.entity_ID)) '(entity ID %s). Coregistration impossible.' % (self.scene_ID, self.entity_ID))
self.spatRef_available = False
return None
# assign spatial refernce scene
self.spatRef_scene = RSF.choose_ref_scene()
self.spatRef_available = True
self.logger.info('Found a suitable reference image for coregistration: scene ID %s (entity ID %s).'
% (self.spatRef_scene.scene_ID, self.spatRef_scene.entity_ID))
def _get_reference_image_params_pgSQL(self): def _get_reference_image_params_pgSQL(self):
# TODO implement earlier version of this function as a backup for SpatialIndexMediator # TODO implement earlier version of this function as a backup for SpatialIndexMediator
......
...@@ -156,16 +156,16 @@ def get_info_from_postgreSQLdb(conn_params, tablename, vals2return, cond_dict=No ...@@ -156,16 +156,16 @@ def get_info_from_postgreSQLdb(conn_params, tablename, vals2return, cond_dict=No
:param tablename: <str> name of the table within the database to be queried :param tablename: <str> name of the table within the database to be queried
:param vals2return: <list or str> a list of strings containing the column titles of the values to be returned :param vals2return: <list or str> a list of strings containing the column titles of the values to be returned
:param cond_dict: <dict> a dictionary containing the query conditions in the form {'column_name':<value>} :param cond_dict: <dict> a dictionary containing the query conditions in the form {'column_name':<value>}
HINT: <value> can also be a list or a tuple of elements to match HINT: <value> can also be a list or a tuple of elements to match, BUT note that the order
of the list items is NOT respected!
:param records2fetch: <int> number of records to be fetched (default=0: fetch unlimited records) :param records2fetch: <int> number of records to be fetched (default=0: fetch unlimited records)
:param timeout: <int> allows to set a custom statement timeout (milliseconds) :param timeout: <int> allows to set a custom statement timeout (milliseconds)
""" """
if not isinstance(vals2return, list): if not isinstance(vals2return, list):
vals2return = [vals2return] vals2return = [vals2return]
assert isinstance(records2fetch, int), \ assert isinstance(records2fetch, int), "get_info_from_postgreSQLdb: Expected an integer for the argument " \
"get_info_from_postgreSQLdb: Expected an integer for the argument 'records2return'. Got %s" % type( "'records2return'. Got %s" % type(records2fetch)
records2fetch)
cond_dict = cond_dict if cond_dict else {} cond_dict = cond_dict if cond_dict else {}
conn_params = "%s options = '-c statement_timeout=%s'" % (conn_params, timeout) conn_params = "%s options = '-c statement_timeout=%s'" % (conn_params, timeout)
connection = psycopg2.connect(conn_params) connection = psycopg2.connect(conn_params)
......
...@@ -802,7 +802,7 @@ class GMS_object(Dataset): ...@@ -802,7 +802,7 @@ class GMS_object(Dataset):
sampleTile = dict(tiles[0]) sampleTile = dict(tiles[0])
target_shape = self.shape_fullArr if len(sampleTile['data'].shape) == 3 else self.shape_fullArr[:2] target_shape = self.shape_fullArr if len(sampleTile['data'].shape) == 3 else self.shape_fullArr[:2]
setattr(self, target_attr, np.empty(target_shape, dtype=sampleTile['data'].dtype)) setattr(self, target_attr, np.empty(target_shape, dtype=sampleTile['data'].dtype))
for tile in tiles: for tile in tiles: # type: dict
rS, rE, cS, cE = tile['row_start'], tile['row_end'], tile['col_start'], tile['col_end'] rS, rE, cS, cE = tile['row_start'], tile['row_end'], tile['col_start'], tile['col_end']
if len(target_shape) == 3: if len(target_shape) == 3:
getattr(self, target_attr)[rS:rE + 1, cS:cE + 1, :] = tile['data'] getattr(self, target_attr)[rS:rE + 1, cS:cE + 1, :] = tile['data']
......
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