Commit 1a5c6148 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

FMASK is now called from atmospheric correction.

parent 6b215624
Pipeline #441 passed with stage
in 5 minutes and 14 seconds
......@@ -45,7 +45,7 @@ class L1A_object(GMS_object):
# TODO docstring
# load default attribute values and methods
super().__init__()
super(L1A_object, self).__init__()
# unpack kwargs
self.proc_level = 'L1A'
......@@ -687,9 +687,7 @@ class L1A_object(GMS_object):
def calc_cloud_mask(self, subset=None):
algorithm = CFG.job.cloud_masking_algorithm
if algorithm not in ['FMASK', 'Classical Bayesian']: # TODO move validation to config
raise ValueError(algorithm)
algorithm = CFG.job.cloud_masking_algorithm[self.satellite]
(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):
# 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]
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. '
'Please make sure that the produced cloud mask does not use potential fill values. '
'Otherwise cloud mask results can be affected.' %outFill)
......
......@@ -32,6 +32,7 @@ from .L1B_P import L1B_object
from ..model.METADATA import get_LayerBandsAssignment
from ..misc.definition_dicts import get_outFillZeroSaturated, proc_chain, get_mask_classdefinition
from ..io.Input_reader import SRF
# from .cloud_masking import Cloud_Mask_Creator # circular dependencies
from sicor.sicor_ac import ac_gms
from sicor.sensors import RSImage
......@@ -681,31 +682,6 @@ class AtmCorr(object):
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):
"""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):
# check if input GMS objects provide a cloud mask
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
if list(set(avail_cloud_masks.values())) == [None]:
inObj2use.mask_clouds = self._calc_mask_clouds(tgt_res)
cm_array = cm_gA[:]
cm_legend = cm_gA.legend
cm_geocoding = cm_gA.gt
return None
if no_avail_CMs:
algorithm = CFG.job.cloud_masking_algorithm[self.inObjs[0].satellite]
# check cloud mask target resolution
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]
if algorithm == 'SICOR':
return None
# get mask array
cm_array = inObj2use.masks[:,:,1] # FIXME hardcoded
else:
# 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
cm_legend = get_mask_classdefinition('mask_clouds')
#{'Clear': 10, 'Thick Clouds': 20, 'Thin Clouds': 30, 'Snow': 40} # FIXME hardcoded
else:
# check if there is a cloud mask with suitable GSD
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
# 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]
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,
mask_legend = cm_legend,
......
......@@ -10,6 +10,7 @@ import tarfile
import zipfile
import itertools
import warnings
import re
# custom
try:
......@@ -21,6 +22,7 @@ import fmask
from ..misc.helper_functions import convert_absPathArchive_to_GDALvsiPath, subcall_with_output
from ..misc.definition_dicts import get_mask_classdefinition
from ..misc.database_tools import get_info_from_postgreSQLdb
from ..config import is_GMSConfig_available
from geoarray import GeoArray
......@@ -117,7 +119,7 @@ class _FMASK_Runner(object):
if self.is_GMSConfig_available:
self.cloud_mask.legend = \
get_mask_classdefinition('mask_clouds')
get_mask_classdefinition('mask_clouds', self.satellite)
else: # use default FMASK legend
warnings.warn('GMS configuration not available. Using default cloud mask legend.')
self.cloud_mask.legend = \
......@@ -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:
self.profiling = profiling
self.benchmark_global = bench_all
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.export_VZA_SZA_SAA_RAA_stats = True # hardcoded
self.export_L1C_obj_dumps = False # hardcoded
......@@ -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):
def mask_to_ENVI_Classification(InObj,maskname):
# type: (object,str) -> (np.ndarray, dict, list, list)
cd = get_mask_classdefinition(maskname)
cd = get_mask_classdefinition(maskname, InObj.satellite)
if cd is None:
InObj.logger.warning("No class definition available for mask '%s'." % maskname)
else:
......
......@@ -75,7 +75,7 @@ def get_GMS_sensorcode(GMS_identifier):
raise
def get_mask_classdefinition(maskname):
def get_mask_classdefinition(maskname, satellite):
if maskname== 'mask_nodata':
return {'No data': 0,
'Data' : 1 }
......@@ -102,7 +102,7 @@ def get_mask_classdefinition(maskname):
'Snow' : 60} # SICOR
}
return legends[CFG.job.cloud_masking_algorithm]
return legends[CFG.job.cloud_masking_algorithm[satellite]]
else:
raise ValueError("'%s' is not a supported mask name." %maskname)
......
......@@ -64,6 +64,7 @@ initArgsDict = {'L1A': (None,), 'L1B': (None,), 'L1C': (None,),
def silentremove(filename):
# type: (str) -> None
"""Remove the given file without raising OSError exceptions, e.g. if the file does not exist."""
try:
os.remove(filename)
except OSError as e:
......@@ -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.
: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)
out, err = proc.communicate()
exitcode = proc.returncode
......@@ -116,6 +118,7 @@ def sorted_nicely(iterable):
""" Sort the given iterable in the way that humans expect.
:param iterable:
"""
convert = lambda text: int(text) if text.isdigit() else text
alphanum_key = lambda key: [ convert(c) for c in re.split('([0-9]+)', key) ]
return sorted(iterable, key = alphanum_key)
......@@ -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 shape_fullArr: <tuple> (rows,columns,bands)
"""
if shape_fullArr is None:
shape_fullArr = json.load(open(path_GMS_file))['shape_fullArr']
rows,cols = shape_fullArr[:2]
......@@ -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()
"""
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)
......@@ -232,6 +237,7 @@ def merge_GMS_tiles_to_GMS_obj(list_GMS_tiles):
@jit
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"""
for arrname in list_arraynames:
samplearray = getattr(list_GMS_tiles[0], arrname)
assert isinstance(samplearray,np.ndarray),\
......@@ -287,6 +293,7 @@ def mp_initializer(globals, globs):
:param globals:
:param globs:
"""
for name, value in globs.items():
try:
value._init(globals, name)
......@@ -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.
:param CornerLonLat: list of XY-coordinate tuples
"""
return str(Polygon(CornerLonLat))
......@@ -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].
:param pGSQL_poly:
"""
if not pGSQL_poly.startswith('POLYGON'):
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)]
......@@ -375,6 +384,7 @@ def get_imageCoords_from_shapelyPoly(shapelyPoly, im_gt):
:param shapelyPoly: <shapely.Polygon>
: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)
coordsArr = get_coordsArr(shapelyPoly)
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):
: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)
"""
rows, cols = arr_shape[:2]
xmin, xmax, ymin, ymax = tgt_bounds
if buffer:
......@@ -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
image bounding box
"""
shpPolyImPrj = reproject_shapelyGeometry(shpPolyLonLat, 4326, im_prj)
imCoordsXY = get_imageCoords_from_shapelyPoly(shpPolyImPrj,im_gt)
bounds = corner_coord_to_minmax(imCoordsXY)
......@@ -437,6 +449,7 @@ def get_UL_LR_from_shapefile_features(path_shp):
:param path_shp: <str> the path of the shapefile
"""
dataSource = ogr.Open(path_shp)
layer = dataSource.GetLayer(0)
ullr_list = []
......@@ -445,11 +458,13 @@ def get_UL_LR_from_shapefile_features(path_shp):
ul = e[0], e[3]
lr = e[1], e[2]
ullr_list.append((ul,lr))
dataSource = layer = None
return ullr_list
def reorder_CornerLonLat(CornerLonLat):
"""Reorders corner coordinate lists from [UL,UR,LL,LR] to clockwise order: [UL,UR,LR,LL]"""
if len(CornerLonLat) > 4:
warnings.warn('Only 4 of the given %s corner coordinates were respected.' %len(CornerLonLat))
return [CornerLonLat[0], CornerLonLat[1], CornerLonLat[3], CornerLonLat[2] ]
......@@ -458,10 +473,12 @@ def reorder_CornerLonLat(CornerLonLat):
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
the database geometry field. """
try:
pgSQL_geom = DB_T.get_info_from_postgreSQLdb(CFG.job.conn_database,'scenes_proc','bounds', {'sceneid':scene_ID})[0][0]
except IndexError:
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,\
'Database error: Received an invalid geometry from the postgreSQL database!'
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):
:param findall: if given, at final a findall
:return: found xml object, None if nothing was found
"""
buf = xml_root.find(str(QName(namespace, branch)))
for br in branches:
buf = buf.find(br)
......@@ -506,6 +524,7 @@ def find_in_xml(xml, *branch):
:param branch: iterate to branches using find
:return: xml object, None if nothing was found
"""
buf = xml
for br in branch:
buf = buf.find(br)
......@@ -519,6 +538,7 @@ def get_values_from_xml(leaf, dtype=np.float):
:param dtype: dtype of returned numpy array
:return: numpy array
"""
return np.array([ele.text.split(" ") for ele in leaf.findall("VALUES")], dtype=dtype)
......
......@@ -13,7 +13,7 @@ import os
from geoarray import GeoArray
from geomultisens import __file__
from geomultisens.algorithms.fmask_runner import FMASK_Runner_Landsat, FMASK_Runner_Sentinel2
from geomultisens.algorithms.cloud_masking import FMASK_Runner_Landsat, FMASK_Runner_Sentinel2
......
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