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

FMASK is now called from atmospheric correction.


Former-commit-id: 1a5c6148
parent d6b6f163
...@@ -45,7 +45,7 @@ class L1A_object(GMS_object): ...@@ -45,7 +45,7 @@ class L1A_object(GMS_object):
# TODO docstring # TODO docstring
# load default attribute values and methods # load default attribute values and methods
super().__init__() super(L1A_object, self).__init__()
# unpack kwargs # unpack kwargs
self.proc_level = 'L1A' self.proc_level = 'L1A'
...@@ -687,9 +687,7 @@ class L1A_object(GMS_object): ...@@ -687,9 +687,7 @@ class L1A_object(GMS_object):
def calc_cloud_mask(self, subset=None): def calc_cloud_mask(self, subset=None):
algorithm = CFG.job.cloud_masking_algorithm algorithm = CFG.job.cloud_masking_algorithm[self.satellite]
if algorithm not in ['FMASK', 'Classical Bayesian']: # TODO move validation to config
raise ValueError(algorithm)
(rS, rE), (cS, cE) = self.arr_pos if self.arr_pos else ((0, self.shape_fullArr[0]), (0, self.shape_fullArr[1])) (rS, rE), (cS, cE) = self.arr_pos if self.arr_pos else ((0, self.shape_fullArr[0]), (0, self.shape_fullArr[1]))
...@@ -820,7 +818,7 @@ class L1A_object(GMS_object): ...@@ -820,7 +818,7 @@ class L1A_object(GMS_object):
# update fill values (fill values written by CLD_obj(self) are not equal to GMS fill values) # update fill values (fill values written by CLD_obj(self) are not equal to GMS fill values)
outFill = get_outFillZeroSaturated(self.mask_clouds.dtype)[0] outFill = get_outFillZeroSaturated(self.mask_clouds.dtype)[0]
if outFill in mask_clouds and get_mask_classdefinition('mask_clouds')['No Data'] != 0: if outFill in mask_clouds and get_mask_classdefinition('mask_clouds', self.satellite)['No Data'] != 0:
warnings.warn('Value %s is assigned to mask_clouds although it already contains this value. ' warnings.warn('Value %s is assigned to mask_clouds although it already contains this value. '
'Please make sure that the produced cloud mask does not use potential fill values. ' 'Please make sure that the produced cloud mask does not use potential fill values. '
'Otherwise cloud mask results can be affected.' %outFill) 'Otherwise cloud mask results can be affected.' %outFill)
......
...@@ -32,6 +32,7 @@ from .L1B_P import L1B_object ...@@ -32,6 +32,7 @@ from .L1B_P import L1B_object
from ..model.METADATA import get_LayerBandsAssignment from ..model.METADATA import get_LayerBandsAssignment
from ..misc.definition_dicts import get_outFillZeroSaturated, proc_chain, get_mask_classdefinition from ..misc.definition_dicts import get_outFillZeroSaturated, proc_chain, get_mask_classdefinition
from ..io.Input_reader import SRF from ..io.Input_reader import SRF
# from .cloud_masking import Cloud_Mask_Creator # circular dependencies
from sicor.sicor_ac import ac_gms from sicor.sicor_ac import ac_gms
from sicor.sensors import RSImage from sicor.sensors import RSImage
...@@ -681,31 +682,6 @@ class AtmCorr(object): ...@@ -681,31 +682,6 @@ class AtmCorr(object):
return SRF(GMS_identifier_fullScene, wvl_unit='nanometers', format_bandnames=True) return SRF(GMS_identifier_fullScene, wvl_unit='nanometers', format_bandnames=True)
def _calc_mask_clouds(self, target_res):
inObj2use = self.inObjs[0]
if re.search('Landsat', inObj2use.satellite, re.I):
from .fmask_runner import FMASK_Runner_Landsat
FMR = FMASK_Runner_Landsat(inObj2use.path_archive, inObj2use.satellite)
cm_gA = FMR.calc_cloudMask()
mask_clouds = cm_gA[:]
elif re.search('Sentinel-2', inObj2use.satellite, re.I):
from .fmask_runner import FMASK_Runner_Sentinel2
FMR = FMASK_Runner_Sentinel2(inObj2use.path_archive, inObj2use.satellite, scene_ID=inObj2use.scene_ID,
target_res=target_res)
cm_gA = FMR.calc_cloudMask()
mask_clouds = cm_gA[:]
else:
self.logger.error(
'Cloud masking is not yet implemented for %s %s...' % (inObj2use.satellite, inObj2use.sensor))
mask_clouds = None
return mask_clouds
def _get_mask_clouds(self): def _get_mask_clouds(self):
"""Returns an instance of S2Mask in case cloud mask is given by input GMS objects. Otherwise None is returned. """Returns an instance of S2Mask in case cloud mask is given by input GMS objects. Otherwise None is returned.
...@@ -716,39 +692,58 @@ class AtmCorr(object): ...@@ -716,39 +692,58 @@ class AtmCorr(object):
# check if input GMS objects provide a cloud mask # check if input GMS objects provide a cloud mask
avail_cloud_masks = {inObj.GMS_identifier['Subsystem']: inObj.mask_clouds for inObj in self.inObjs} avail_cloud_masks = {inObj.GMS_identifier['Subsystem']: inObj.mask_clouds for inObj in self.inObjs}
no_avail_CMs = list(set(avail_cloud_masks.values())) == [None]
# compute cloud mask if not already provided # compute cloud mask if not already provided
if list(set(avail_cloud_masks.values())) == [None]: if no_avail_CMs:
inObj2use.mask_clouds = self._calc_mask_clouds(tgt_res) algorithm = CFG.job.cloud_masking_algorithm[self.inObjs[0].satellite]
cm_array = cm_gA[:]
cm_legend = cm_gA.legend
cm_geocoding = cm_gA.gt
return None
# check cloud mask target resolution if algorithm == 'SICOR':
inObjs2use = [obj for obj in self.inObjs if obj.mask_clouds is not None and obj.mask_clouds.xgsd==tgt_res] return None
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 else:
cm_array = inObj2use.masks[:,:,1] # FIXME hardcoded # FMASK or Classical Bayesian
try:
from .cloud_masking import Cloud_Mask_Creator
CMC = Cloud_Mask_Creator(self.inObjs[0], algorithm=algorithm)
CMC.calc_cloud_mask()
cm_geoarray = CMC.cloud_mask_geoarray
cm_array = CMC.cloud_mask_array
cm_legend = CMC.cloud_mask_legend
except Exception as err:
print('\nAn error occurred during FMASK cloud masking.')
print("Print traceback in case you care:")
print(traceback.format_exc())
return None
# get legend else:
cm_legend = get_mask_classdefinition('mask_clouds') # check if there is a cloud mask with suitable GSD
#{'Clear': 10, 'Thick Clouds': 20, 'Thin Clouds': 30, 'Snow': 40} # FIXME hardcoded 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 (geo)array
cm_geoarray = inObj2use.mask_clouds
cm_array = inObj2use.mask_clouds[:]
# get legend
cm_legend = get_mask_classdefinition('mask_clouds', inObj2use.satellite)
#{'Clear': 10, 'Thick Clouds': 20, 'Thin Clouds': 30, 'Snow': 40} # FIXME hardcoded
# validate that xGSD equals yGSD
if cm_geoarray.xgsd != cm_geoarray.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.')
# get geocoding # get geocoding
# validate that xGSD equals yGSD
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] cm_geocoding = self.metadata["spatial_samplings"][tgt_res]
self.options['cld_mask']['nodata_value_mask'] = inObj2use.mask_clouds.nodata # get nodata value
self.options['cld_mask']['nodata_value_mask'] = cm_geoarray.nodata
return S2Mask(mask_array = cm_array, return S2Mask(mask_array = cm_array,
mask_legend = cm_legend, mask_legend = cm_legend,
......
...@@ -10,6 +10,7 @@ import tarfile ...@@ -10,6 +10,7 @@ import tarfile
import zipfile import zipfile
import itertools import itertools
import warnings import warnings
import re
# custom # custom
try: try:
...@@ -21,6 +22,7 @@ import fmask ...@@ -21,6 +22,7 @@ import fmask
from ..misc.helper_functions import convert_absPathArchive_to_GDALvsiPath, subcall_with_output from ..misc.helper_functions import convert_absPathArchive_to_GDALvsiPath, subcall_with_output
from ..misc.definition_dicts import get_mask_classdefinition from ..misc.definition_dicts import get_mask_classdefinition
from ..misc.database_tools import get_info_from_postgreSQLdb from ..misc.database_tools import get_info_from_postgreSQLdb
from ..config import is_GMSConfig_available
from geoarray import GeoArray from geoarray import GeoArray
...@@ -117,7 +119,7 @@ class _FMASK_Runner(object): ...@@ -117,7 +119,7 @@ class _FMASK_Runner(object):
if self.is_GMSConfig_available: if self.is_GMSConfig_available:
self.cloud_mask.legend = \ self.cloud_mask.legend = \
get_mask_classdefinition('mask_clouds') get_mask_classdefinition('mask_clouds', self.satellite)
else: # use default FMASK legend else: # use default FMASK legend
warnings.warn('GMS configuration not available. Using default cloud mask legend.') warnings.warn('GMS configuration not available. Using default cloud mask legend.')
self.cloud_mask.legend = \ self.cloud_mask.legend = \
...@@ -361,3 +363,65 @@ class FMASK_Runner_Sentinel2(_FMASK_Runner): ...@@ -361,3 +363,65 @@ class FMASK_Runner_Sentinel2(_FMASK_Runner):
class Cloud_Mask_Creator(object):
def __init__(self, GMS_object, algorithm, target_res=None):
"""A class for creating cloud masks.
:param GMS_object: instance or subclass instance of model.gms_object.GMS_object
:param algorithm: 'FMASK' or 'Classical Bayesian
:param target_res: target resolution of the computed cloud mask (if not given, the appropriate resolution
needed for atmospheric correction is chosen)
"""
self.GMS_obj = GMS_object
self.algorithm = algorithm
self.tgt_res = target_res if target_res else self.GMS_obj.ac_options['cld_mask']['target_resolution']
self.cloud_mask_geoarray = None
self.cloud_mask_array = None
self.cloud_mask_legend = None
self.success = None
if not self.GMS_obj.satellite in ['Landsat-4', 'Landsat-5', 'Landsat-7', 'Landsat-8',
'Sentinel-2A', 'Sentinel-2B']:
self.GMS_obj.error(
'Cloud masking is not yet implemented for %s %s...' % (self.GMS_obj.satellite, self.GMS_obj.sensor))
self.success = False
if algorithm not in ['FMASK', 'Classical Bayesian']: # TODO move validation to config
raise ValueError(algorithm)
def calc_cloud_mask(self):
if self.success is False:
return None, None, None
self.GMS_obj.logger.info("Calculating cloud mask based on '%s' algorithm..." % self.algorithm)
if self.algorithm is 'FMASK':
if re.search('Landsat', self.GMS_obj.satellite, re.I):
FMR = FMASK_Runner_Landsat(self.GMS_obj.path_archive, self.GMS_obj.satellite)
else:
FMR = FMASK_Runner_Sentinel2(self.GMS_obj.path_archive, self.GMS_obj.satellite,
scene_ID = self.GMS_obj.scene_ID,
target_res = self.tgt_res)
self.cloud_mask_geoarray = FMR.calc_cloudMask()
self.cloud_mask_array = self.cloud_mask_geoarray[:]
self.cloud_mask_legend = self.cloud_mask_geoarray.legend
else:
raise NotImplementedError("The cloud masking algorithm '%s' is not yet implemented." %self.algorithm)
self.validate_cloud_mask()
return self.cloud_mask_geoarray, self.cloud_mask_array, self.cloud_mask_legend
def validate_cloud_mask(self):
assert self.cloud_mask_geoarray.xgsd == self.cloud_mask_geoarray.ygsd == self.tgt_res
self.success = True
...@@ -134,7 +134,13 @@ class Job: ...@@ -134,7 +134,13 @@ class Job:
self.profiling = profiling self.profiling = profiling
self.benchmark_global = bench_all self.benchmark_global = bench_all
self.bench_CLD_class = bench_cloudMask self.bench_CLD_class = bench_cloudMask
self.cloud_masking_algorithm = 'FMASK' # 'FMASK', 'Classical Bayesian', 'SICOR' self.cloud_masking_algorithm = {'Landsat-4': 'FMASK',
'Landsat-5': 'FMASK',
'Landsat-7': 'FMASK',
'Landsat-8': 'FMASK',
'Sentinel-2A': 'FMASK',
'Sentinel-2B': 'FMASK',
} # 'FMASK', 'Classical Bayesian', 'SICOR'
self.SZA_SAA_calculation_accurracy = 'coarse' # hardcoded self.SZA_SAA_calculation_accurracy = 'coarse' # hardcoded
self.export_VZA_SZA_SAA_RAA_stats = True # hardcoded self.export_VZA_SZA_SAA_RAA_stats = True # hardcoded
self.export_L1C_obj_dumps = False # hardcoded self.export_L1C_obj_dumps = False # hardcoded
...@@ -510,3 +516,9 @@ class Usecase: ...@@ -510,3 +516,9 @@ class Usecase:
def is_GMSConfig_available():
try:
if GMS_config.job is not None:
return True
except (EnvironmentError, OSError):
return False
...@@ -290,7 +290,7 @@ def reorder_ENVI_header(path_hdr,tgt_keyOrder): ...@@ -290,7 +290,7 @@ def reorder_ENVI_header(path_hdr,tgt_keyOrder):
def mask_to_ENVI_Classification(InObj,maskname): def mask_to_ENVI_Classification(InObj,maskname):
# type: (object,str) -> (np.ndarray, dict, list, list) # type: (object,str) -> (np.ndarray, dict, list, list)
cd = get_mask_classdefinition(maskname) cd = get_mask_classdefinition(maskname, InObj.satellite)
if cd is None: if cd is None:
InObj.logger.warning("No class definition available for mask '%s'." % maskname) InObj.logger.warning("No class definition available for mask '%s'." % maskname)
else: else:
......
...@@ -75,7 +75,7 @@ def get_GMS_sensorcode(GMS_identifier): ...@@ -75,7 +75,7 @@ def get_GMS_sensorcode(GMS_identifier):
raise raise
def get_mask_classdefinition(maskname): def get_mask_classdefinition(maskname, satellite):
if maskname== 'mask_nodata': if maskname== 'mask_nodata':
return {'No data': 0, return {'No data': 0,
'Data' : 1 } 'Data' : 1 }
...@@ -102,7 +102,7 @@ def get_mask_classdefinition(maskname): ...@@ -102,7 +102,7 @@ def get_mask_classdefinition(maskname):
'Snow' : 60} # SICOR 'Snow' : 60} # SICOR
} }
return legends[CFG.job.cloud_masking_algorithm] return legends[CFG.job.cloud_masking_algorithm[satellite]]
else: else:
raise ValueError("'%s' is not a supported mask name." %maskname) raise ValueError("'%s' is not a supported mask name." %maskname)
......
...@@ -64,6 +64,7 @@ initArgsDict = {'L1A': (None,), 'L1B': (None,), 'L1C': (None,), ...@@ -64,6 +64,7 @@ initArgsDict = {'L1A': (None,), 'L1B': (None,), 'L1C': (None,),
def silentremove(filename): def silentremove(filename):
# type: (str) -> None # type: (str) -> None
"""Remove the given file without raising OSError exceptions, e.g. if the file does not exist.""" """Remove the given file without raising OSError exceptions, e.g. if the file does not exist."""
try: try:
os.remove(filename) os.remove(filename)
except OSError as e: except OSError as e:
...@@ -105,6 +106,7 @@ def subcall_with_output(cmd, no_stdout=False, no_stderr=False): ...@@ -105,6 +106,7 @@ def subcall_with_output(cmd, no_stdout=False, no_stderr=False):
"""Execute external command and get its stdout, exitcode and stderr. """Execute external command and get its stdout, exitcode and stderr.
:param cmd: a normal shell command including parameters :param cmd: a normal shell command including parameters
""" """
proc = Popen(shlex.split(cmd), stdout=None if no_stdout else PIPE, stderr=None if no_stderr else PIPE) proc = Popen(shlex.split(cmd), stdout=None if no_stdout else PIPE, stderr=None if no_stderr else PIPE)
out, err = proc.communicate() out, err = proc.communicate()
exitcode = proc.returncode exitcode = proc.returncode
...@@ -116,6 +118,7 @@ def sorted_nicely(iterable): ...@@ -116,6 +118,7 @@ def sorted_nicely(iterable):
""" Sort the given iterable in the way that humans expect. """ Sort the given iterable in the way that humans expect.
:param iterable: :param iterable:
""" """
convert = lambda text: int(text) if text.isdigit() else text convert = lambda text: int(text) if text.isdigit() else text
alphanum_key = lambda key: [ convert(c) for c in re.split('([0-9]+)', key) ] alphanum_key = lambda key: [ convert(c) for c in re.split('([0-9]+)', key) ]
return sorted(iterable, key = alphanum_key) return sorted(iterable, key = alphanum_key)
...@@ -146,6 +149,7 @@ def get_image_tileborders(target_tileShape, target_tileSize, path_GMS_file=None, ...@@ -146,6 +149,7 @@ def get_image_tileborders(target_tileShape, target_tileSize, path_GMS_file=None,
:param path_GMS_file: <str> path to the *.gms file. Only needed if shape_fullArr is not given :param path_GMS_file: <str> path to the *.gms file. Only needed if shape_fullArr is not given
:param shape_fullArr: <tuple> (rows,columns,bands) :param shape_fullArr: <tuple> (rows,columns,bands)
""" """
if shape_fullArr is None: if shape_fullArr is None:
shape_fullArr = json.load(open(path_GMS_file))['shape_fullArr'] shape_fullArr = json.load(open(path_GMS_file))['shape_fullArr']
rows,cols = shape_fullArr[:2] rows,cols = shape_fullArr[:2]
...@@ -211,6 +215,7 @@ def merge_GMS_tiles_to_GMS_obj(list_GMS_tiles): ...@@ -211,6 +215,7 @@ def merge_GMS_tiles_to_GMS_obj(list_GMS_tiles):
:param list_GMS_tiles: <list> of GMS objects that have been created by cut_GMS_obj_into_blocks() :param list_GMS_tiles: <list> of GMS objects that have been created by cut_GMS_obj_into_blocks()
""" """
warnings.warn("'merge_GMS_tiles_to_GMS_obj' is deprecated. Use GMS_object.from_tiles() instead!", DeprecationWarning) warnings.warn("'merge_GMS_tiles_to_GMS_obj' is deprecated. Use GMS_object.from_tiles() instead!", DeprecationWarning)
if 'IMapUnorderedIterator' in str(type(list_GMS_tiles)): list_GMS_tiles = list(list_GMS_tiles) if 'IMapUnorderedIterator' in str(type(list_GMS_tiles)): list_GMS_tiles = list(list_GMS_tiles)
...@@ -232,6 +237,7 @@ def merge_GMS_tiles_to_GMS_obj(list_GMS_tiles): ...@@ -232,6 +237,7 @@ def merge_GMS_tiles_to_GMS_obj(list_GMS_tiles):
@jit @jit
def _numba_array_merger(GMS_obj, list_arraynames, list_GMS_tiles): def _numba_array_merger(GMS_obj, list_arraynames, list_GMS_tiles):
"""private function, e.g. called by merge_GMS_tiles_to_GMS_obj() in order to fasten array merging""" """private function, e.g. called by merge_GMS_tiles_to_GMS_obj() in order to fasten array merging"""
for arrname in list_arraynames: for arrname in list_arraynames:
samplearray = getattr(list_GMS_tiles[0], arrname) samplearray = getattr(list_GMS_tiles[0], arrname)
assert isinstance(samplearray,np.ndarray),\ assert isinstance(samplearray,np.ndarray),\
...@@ -287,6 +293,7 @@ def mp_initializer(globals, globs): ...@@ -287,6 +293,7 @@ def mp_initializer(globals, globs):
:param globals: :param globals:
:param globs: :param globs:
""" """
for name, value in globs.items(): for name, value in globs.items():
try: try:
value._init(globals, name) value._init(globals, name)
...@@ -331,6 +338,7 @@ def cornerLonLat_to_postgreSQL_poly(CornerLonLat): ...@@ -331,6 +338,7 @@ def cornerLonLat_to_postgreSQL_poly(CornerLonLat):
"""Converts a coordinate list [UL_LonLat, UR_LonLat, LL_LonLat, LR_LonLat] to a postgreSQL polygon. """Converts a coordinate list [UL_LonLat, UR_LonLat, LL_LonLat, LR_LonLat] to a postgreSQL polygon.
:param CornerLonLat: list of XY-coordinate tuples :param CornerLonLat: list of XY-coordinate tuples
""" """
return str(Polygon(CornerLonLat)) return str(Polygon(CornerLonLat))
...@@ -339,6 +347,7 @@ def postgreSQL_poly_to_cornerLonLat(pGSQL_poly): ...@@ -339,6 +347,7 @@ def postgreSQL_poly_to_cornerLonLat(pGSQL_poly):
"""Converts a postgreSQL polygon to a coordinate list [UL_LonLat, UR_LonLat, LL_LonLat, LR_LonLat]. """Converts a postgreSQL polygon to a coordinate list [UL_LonLat, UR_LonLat, LL_LonLat, LR_LonLat].
:param pGSQL_poly: :param pGSQL_poly:
""" """
if not pGSQL_poly.startswith('POLYGON'): if not pGSQL_poly.startswith('POLYGON'):
raise ValueError("'pGSQL_poly' has to start with 'POLYGON...'. Got %s" %pGSQL_poly) raise ValueError("'pGSQL_poly' has to start with 'POLYGON...'. Got %s" %pGSQL_poly)
fl = [float(i) for i in re.findall(r"[-+]?\d*\.\d+|\d+",pGSQL_poly)] fl = [float(i) for i in re.findall(r"[-+]?\d*\.\d+|\d+",pGSQL_poly)]
...@@ -375,6 +384,7 @@ def get_imageCoords_from_shapelyPoly(shapelyPoly, im_gt): ...@@ -375,6 +384,7 @@ def get_imageCoords_from_shapelyPoly(shapelyPoly, im_gt):
:param shapelyPoly: <shapely.Polygon> :param shapelyPoly: <shapely.Polygon>
:param im_gt: <list> the GDAL geotransform of the target image :param im_gt: <list> the GDAL geotransform of the target image
""" """
get_coordsArr = lambda shpPoly: np.swapaxes(np.array(shpPoly.exterior.coords.xy), 0, 1) get_coordsArr = lambda shpPoly: np.swapaxes(np.array(shpPoly.exterior.coords.xy), 0, 1)
coordsArr = get_coordsArr(shapelyPoly) coordsArr = get_coordsArr(shapelyPoly)
imCoordsXY = [mapXY2imXY((X, Y), im_gt) for X, Y in coordsArr.tolist()] imCoordsXY = [mapXY2imXY((X, Y), im_gt) for X, Y in coordsArr.tolist()]
...@@ -392,6 +402,7 @@ def get_valid_arrSubsetBounds(arr_shape, tgt_bounds, buffer=0): ...@@ -392,6 +402,7 @@ def get_valid_arrSubsetBounds(arr_shape, tgt_bounds, buffer=0):
:param tgt_bounds: <tuple of floats> the target image coordinates in the form (xmin, xmax, ymin, ymax) :param tgt_bounds: <tuple of floats> the target image coordinates in the form (xmin, xmax, ymin, ymax)
:param buffer: <float> an optional buffer size (image pixel units) :param buffer: <float> an optional buffer size (image pixel units)
""" """
rows, cols = arr_shape[:2] rows, cols = arr_shape[:2]
xmin, xmax, ymin, ymax = tgt_bounds xmin, xmax, ymin, ymax = tgt_bounds
if buffer: if buffer:
...@@ -420,6 +431,7 @@ def get_arrSubsetBounds_from_shapelyPolyLonLat(arr_shape, shpPolyLonLat, im_gt, ...@@ -420,6 +431,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 :param ensure_valid_coords: <bool> whether to ensure that the returned values are all inside the original
image bounding box image bounding box
""" """
shpPolyImPrj = reproject_shapelyGeometry(shpPolyLonLat, 4326, im_prj) shpPolyImPrj = reproject_shapelyGeometry(shpPolyLonLat, 4326, im_prj)
imCoordsXY = get_imageCoords_from_shapelyPoly(shpPolyImPrj,im_gt) imCoordsXY = get_imageCoords_from_shapelyPoly(shpPolyImPrj,im_gt)
bounds = corner_coord_to_minmax(imCoordsXY) bounds = corner_coord_to_minmax(imCoordsXY)
...@@ -437,6 +449,7 @@ def get_UL_LR_from_shapefile_features(path_shp): ...@@ -437,6 +449,7 @@ def get_UL_LR_from_shapefile_features(path_shp):
:param path_shp: <str> the path of the shapefile :param path_shp: <str> the path of the shapefile
""" """
dataSource = ogr.Open(path_shp) dataSource = ogr.Open(path_shp)
layer = dataSource.GetLayer(0) layer = dataSource.GetLayer(0)
ullr_list = [] ullr_list = []
...@@ -445,11 +458,13 @@ def get_UL_LR_from_shapefile_features(path_shp): ...@@ -445,11 +458,13 @@ def get_UL_LR_from_shapefile_features(path_shp):
ul = e[0], e[3] ul = e[0], e[3]
lr = e[1], e[2] lr = e[1], e[2]
ullr_list.append((ul,lr)) ullr_list.append((ul,lr))
dataSource = layer = None
return ullr_list return ullr_list
def reorder_CornerLonLat(CornerLonLat): def reorder_CornerLonLat(CornerLonLat):
"""Reorders corner coordinate lists from [UL,UR,LL,LR] to clockwise order: [UL,UR,LR,LL]""" """Reorders corner coordinate lists from [UL,UR,LL,LR] to clockwise order: [UL,UR,LR,LL]"""
if len(CornerLonLat) > 4: if len(CornerLonLat) > 4:
warnings.warn('Only 4 of the given %s corner coordinates were respected.' %len(CornerLonLat)) warnings.warn('Only 4 of the given %s corner coordinates were respected.' %len(CornerLonLat))
return [CornerLonLat[0], CornerLonLat[1], CornerLonLat[3], CornerLonLat[2] ] return [CornerLonLat[0], CornerLonLat[1], CornerLonLat[3], CornerLonLat[2] ]
...@@ -458,10 +473,12 @@ def reorder_CornerLonLat(CornerLonLat): ...@@ -458,10 +473,12 @@ def reorder_CornerLonLat(CornerLonLat):
def sceneID_to_trueDataCornerLonLat(scene_ID): def sceneID_to_trueDataCornerLonLat(scene_ID):
"""Returns a list of corner coordinates ordered like (UL,UR,LL,LR) corresponding to the given scene_ID by querying """Returns a list of corner coordinates ordered like (UL,UR,LL,LR) corresponding to the given scene_ID by querying
the database geometry field. """ the database geometry field. """
try: try:
pgSQL_geom = DB_T.get_info_from_postgreSQLdb(CFG.job.conn_database,'scenes_proc','bounds', {'sceneid':scene_ID})[0][0] pgSQL_geom = DB_T.get_info_from_postgreSQLdb(CFG.job.conn_database,'scenes_proc','bounds', {'sceneid':scene_ID})[0][0]
except IndexError: except IndexError:
pgSQL_geom = DB_T.get_info_from_postgreSQLdb(CFG.job.conn_database,'scenes','bounds', {'id':scene_ID})[0][0] pgSQL_geom = DB_T.get_info_from_postgreSQLdb(CFG.job.conn_database,'scenes','bounds', {'id':scene_ID})[0][0]
assert shapely.wkb.loads(pgSQL_geom,hex=True).is_valid,\ assert shapely.wkb.loads(pgSQL_geom,hex=True).is_valid,\
'Database error: Received an invalid geometry from the postgreSQL database!' 'Database error: Received an invalid geometry from the postgreSQL database!'
return postgreSQL_poly_to_cornerLonLat(postgreSQL_geometry_to_postgreSQL_poly(pgSQL_geom)) return postgreSQL_poly_to_cornerLonLat(postgreSQL_geometry_to_postgreSQL_poly(pgSQL_geom))
...@@ -491,6 +508,7 @@ def find_in_xml_root(namespace, xml_root, branch, *branches, findall=None): ...@@ -491,6 +508,7 @@ def find_in_xml_root(namespace, xml_root, branch, *branches, findall=None):
:param findall: if given, at final a findall :param findall: if given, at final a findall
:return: found xml object, None if nothing was found :return: found xml object, None if nothing was found
""" """
buf = xml_root.find(str(QName(namespace, branch))) buf = xml_root.find(str(QName(namespace, branch)))
for br in branches: for br in branches:
buf = buf.find(br) buf = buf.find(br)
...@@ -506,6 +524,7 @@ def find_in_xml(xml, *branch): ...@@ -506,6 +524,7 @@ def find_in_xml(xml, *branch):
:param branch: iterate to branches using find :param branch: iterate to branches using find
:return: xml object, None if nothing was found :return: xml object, None if nothing was found
""" """
buf = xml buf = xml
for br in branch: for br in branch:
buf = buf.find(br) buf = buf.find(br)
...@@ -519,6 +538,7 @@ def get_values_from_xml(leaf, dtype=np.float): ...@@ -519,6 +538,7 @@ def get_values_from_xml(leaf, dtype=np.float):
:param dtype: dtype of returned numpy array :param dtype: dtype of returned numpy array
:return: numpy array :return: numpy array
""" """