Commit 0fdfaaed authored by Daniel Scheffler's avatar Daniel Scheffler Committed by Mathias Peters
Browse files

Bugfix for exception when merging multiple sensor subsystems

algorithms.GEOPROCESSING:
- added get_common_extent()
algorithms.L1B_P,L1B_object:
- implemented clipextent and clipextent_prj keywords
processing.pipeline.L2A_map:
- updated correct_spatial_shifts call
- updated __version__
Former-commit-id: 3aec298d
Former-commit-id: c860ecc3
parent ce2e9a10
......@@ -15,7 +15,7 @@ from . import config
from .processing.process_controller import process_controller
__version__ = '20170125.01'
__version__ = '20170125.02'
__author__ = 'Daniel Scheffler'
__all__ = ['algorithms',
'io',
......
......@@ -47,7 +47,8 @@ except ImportError:
import pyproj
from pyorbital import astronomy
import ephem
from shapely.geometry import MultiPoint
from shapely.geometry import MultiPoint, Polygon
from shapely.ops import cascaded_union
from py_tools_ds.ptds import GeoArray
from py_tools_ds.ptds.geo.coord_grid import snap_bounds_to_pixGrid
......@@ -2751,6 +2752,27 @@ def get_point_bounds(points):
return min_lon, min_lat, max_lon, max_lat
def get_common_extent(list_extentcoords, alg='outer', return_box=True):
"""Returns the common extent coordinates of all given coordinate extents.
NOTE: this function asserts that all input coordinates belong to the same projection!
:param list_extentcoords: a list of coordinate sets, e.g. [ [ULxy1,URxy1,LRxy1,LLxy1], ULxy2,URxy2,LRxy2,LLxy2],]
:param alg: 'outer': return the outer Polygon extent
'inner': return the inner Polygon extent
:param return_box: whether to return a coordinate box/envelope of the output Polygon
:return: [ULxy,URxy,LRxy,LLxy]
"""
if alg!='outer':
raise NotImplementedError # TODO apply shapely intersection
allPolys = [Polygon(extentcoords).envelope for extentcoords in list_extentcoords]
common_extent = cascaded_union(allPolys)
get_coordsXY = lambda shpPoly: np.swapaxes(np.array(shpPoly.exterior.coords.xy), 0, 1).tolist()
return get_coordsXY(common_extent) if not return_box else get_coordsXY(common_extent.envelope)
def zoom_2Darray_to_shapeFullArr(arr2D,shapeFullArr,meshwidth=1,subset=None,method='linear'):
assert method in ['linear','nearest']
rS,rE,cS,cE = list(get_subsetProps_from_subsetArg(shapeFullArr,subset).values())[3:7]
......
......@@ -698,16 +698,36 @@ class L1B_object(L1A_object):
self.coreg_info.update({'success' : True if not self.coreg_needed else False}) # False means spatRef not available
def correct_spatial_shifts(self, cliptoextent=True, v=False):
"""Corrects the spatial shifts calculated by self.coregister_spatially()."""
def correct_spatial_shifts(self, cliptoextent=True, clipextent=None, clipextent_prj=None, v=False):
# type: (bool, list, any, bool) -> None
"""Corrects the spatial shifts calculated by self.coregister_spatially().
:param cliptoextent: whether to clip the output to the given extent
:param clipextent: list of XY-coordinate tuples giving the target extent (if not given and cliptoextent is
True, the 'trueDataCornerLonLat' attribute of the GMS object is used
:param clipextent_prj: WKT projection string or EPSG code of the projection for the coordinates in clipextent
:param v:
:return:
"""
cliptoextent = cliptoextent if not clipextent else True # cliptoextent is automatically True if an extent is given
if cliptoextent or self.resamp_needed or (self.coreg_needed and self.coreg_info['success']):
## get target bounds # TODO implement boxObj call instead here
trueDataCornerUTM = [transform_any_prj(EPSG2WKT(4326), self.meta_odict['coordinate system string'],x,y)
for x, y in self.trueDataCornerLonLat]
xmin, xmax, ymin, ymax = corner_coord_to_minmax(trueDataCornerUTM)
mapBounds = box(xmin, ymin, xmax, ymax).bounds
if not clipextent:
trueDataCornerUTM = [transform_any_prj(EPSG2WKT(4326), self.MetaObj.projection, x,y)
for x, y in self.trueDataCornerLonLat]
xmin, xmax, ymin, ymax = corner_coord_to_minmax(trueDataCornerUTM)
mapBounds = box(xmin, ymin, xmax, ymax).bounds
else:
assert clipextent_prj, \
"'clipextent_prj' must be given together with 'clipextent'. Received only 'clipextent'."
clipextent_UTM = clipextent if prj_equal(self.MetaObj.projection, clipextent_prj) else \
[transform_any_prj(clipextent_prj, self.MetaObj.projection, x, y) for x, y in clipextent]
xmin, xmax, ymin, ymax = corner_coord_to_minmax(clipextent_UTM)
mapBounds = box(xmin, ymin, xmax, ymax).bounds
# correct shifts and clip to extent
# ensure self.masks exists (does not exist in Flink mode because in that case self.fill_from_disk() is skipped)
......
......@@ -98,7 +98,7 @@ class Job:
self.exec__L1AP = [1, 1, 1]
self.exec__L1BP = [1, 1, 1]
self.exec__L1CP = [1, 1, 1]
self.exec__L2AP = [1, 1, 1]
self.exec__L2AP = [1, 1, 0]
self.exec__L2BP = [1, 1, 0]
self.exec__L2CP = [1, 1, 0]
......
......@@ -9,6 +9,7 @@ from ..algorithms import L1C_P # Level 1C Processor
from ..algorithms import L2A_P # Level 2A Processor
from ..algorithms import L2B_P # Level 2B Processor
from ..algorithms import L2C_P # Level 2C Processor
from ..algorithms.GEOPROCESSING import get_common_extent
@EXC_H.log_uncaught_exceptions
......@@ -165,7 +166,8 @@ def L2A_map(L1C_objs, block_size=None):
L2A_objs = [L2A_P.L2A_object(L1C_obj) for L1C_obj in L1C_objs]
# correct geometric displacements and resample to target grid
[L2A_obj.correct_spatial_shifts() for L2A_obj in L2A_objs]
common_extent = get_common_extent([obj.trueDataCornerLonLat for obj in L2A_objs]) if len(L2A_objs) > 1 else None
[L2A_obj.correct_spatial_shifts(clipextent=common_extent, clipextent_prj=4326) for L2A_obj in L2A_objs]
# merge multiple subsystems belonging to the same scene ID to a single GMS object
L2A_obj = L2A_P.L2A_object().from_sensor_subsystems(L2A_objs) if len(L2A_objs) > 1 else L2A_objs[0]
......
......@@ -668,7 +668,7 @@ def get_job_summary(list_GMS_objects):
all_sat, all_sen = zip(*[i.split('__') for i in (np.unique(DJS['satellite'] + '__' + DJS['sensor']))])
QJS['satellite'] = all_sat
QJS['sensor'] = all_sen
# count objects of with the same satellite/sensor/sceneid combination
# count objects with the same satellite/sensor/sceneid combination
QJS['count'] = [len(DJS[(DJS['satellite'] == sat) & (DJS['sensor'] == sen)]['scene_ID'].unique()) \
for sat, sen in zip(all_sat, all_sen)]
QJS['proc_successfully'] = [len(DJS[(DJS['satellite'] == sat) &
......
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