Commit 3ba2aa8b authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Fixed issue #76 (Cloud mask within *.masks.bsq contains no data values at non-clear positions).

parent e08e16e9
......@@ -964,11 +964,12 @@ class AtmCorr(object):
# apply the original nodata mask (indicating background values)
surf_refl[np.array(inObj.mask_nodata).astype(np.int8) == 0] = oF_refl
# set AC NaNs to GMS outFill
# NOTE: SICOR result has NaNs at no data positions AND non-clear positions
if self.results.bad_data_value is np.nan:
surf_refl[np.isnan(surf_refl)] = oF_refl
else:
surf_refl[
surf_refl == self.results.bad_data_value] = oF_refl # FIXME meaningful to set AC nans to -9999?
surf_refl[surf_refl == self.results.bad_data_value] = oF_refl
# use inObj.arr setter to generate a GeoArray
inObj.arr = surf_refl.astype(inObj.arr.dtype) # -> int16 (also converts NaNs to 0 if needed
......
......@@ -861,19 +861,19 @@ class GMS_object(Dataset):
# self.mask_nodata and self.mask_clouds will already be set here -> so just recreate it from there
GMS_obj_merged.masks = None
# avoid unequal nodata edges between indiviual layers (resampling artifacts etc.) #
###################################################################################
# avoid unequal nodata edges between individual layers (resampling artifacts etc.) #
####################################################################################
# get pixel values of areas that have not been atmospherically corrected (non-clear)
nonclear_labels = [lbl for lbl in ["Clear", "Snow", "Water", "Shadow", "Cirrus", "Cloud"]
if lbl not in CFG.ac_clear_area_labels]
cloud_mask_legend = DEF_D.get_mask_classdefinition('mask_clouds', GMS_obj_merged.satellite)
nonclear_pixVals = [cloud_mask_legend[lbl] for lbl in nonclear_labels]
# apply cloud mask to image data and all products derived from image data
# (only if image data represents BOA-Ref and cloud areas are not to be filled with TOA-Ref)
if re.search('BOA_Reflectance', GMS_obj_merged.MetaObj.PhysUnit, re.I) and not CFG.ac_fillnonclear_areas:
# get pixel values of areas that have not been atmospherically corrected (non-clear)
nonclear_labels = [lbl for lbl in ["Clear", "Snow", "Water", "Shadow", "Cirrus", "Cloud"]
if lbl not in CFG.ac_clear_area_labels]
cloud_mask_legend = DEF_D.get_mask_classdefinition('mask_clouds', GMS_obj_merged.satellite)
nonclear_pixVals = [cloud_mask_legend[lbl] for lbl in nonclear_labels]
# fill non-clear areas with no data values (for all bands)
for pixVal in nonclear_pixVals:
mask_nonclear = GMS_obj_merged.mask_clouds[:] == pixVal
......@@ -883,15 +883,24 @@ class GMS_object(Dataset):
GMS_obj_merged.ac_errors[mask_nonclear] = \
DEF_D.get_outFillZeroSaturated(GMS_obj_merged.ac_errors.dtype)[0]
# update no data mask
# update no data mask (adds all pixels to no data where at least one band has no data)
GMS_obj_merged.calc_mask_nodata(overwrite=True)
# apply updated nodata mask to array data
for attrname in ['arr', 'ac_errors', 'dem', 'mask_clouds', 'mask_clouds_confidence']:
# apply updated nodata mask to array data (only to no data area, not to non-clear area)
# => this sets all pixels to no data where at least one band has no data
# NOTE: This may cause the cloud mask not exactly match the holes in self.arr so that there may be small
# gaps between the cloud edge and where the image data begins. This is due to resampling, e.g. from
# Sentinel-2 60m grid to Landsat 30 m grid. Nodata mask marks those areas as no data although there is
# no cloud at exactly that pixel. This is ok!
area2replace = np.all([GMS_obj_merged.mask_nodata[:] == 0] +
[GMS_obj_merged.mask_clouds[:] != pixVal for pixVal in nonclear_pixVals], axis=0) \
if GMS_obj_merged.mask_clouds is not None else GMS_obj_merged.mask_nodata[:] == 0
for attrname in ['arr', 'ac_errors']:
attr_val = getattr(GMS_obj_merged, attrname)
if attr_val is not None:
attr_val[GMS_obj_merged.mask_nodata[:] == 0] = DEF_D.get_outFillZeroSaturated(attr_val.dtype)[0]
attr_val[area2replace] = DEF_D.get_outFillZeroSaturated(attr_val.dtype)[0]
setattr(GMS_obj_merged, attrname, attr_val)
# recreate self.masks
......@@ -1074,8 +1083,10 @@ class GMS_object(Dataset):
# set self.masks_meta
nodataVal = DEF_D.get_outFillZeroSaturated(self.masks.dtype)[0]
self.masks_meta = {'map info': self.MetaObj.map_info, 'coordinate system string': self.MetaObj.projection,
'bands': len(arrays2combine), 'band names': arrays2combine,
self.masks_meta = {'map info': self.MetaObj.map_info,
'coordinate system string': self.MetaObj.projection,
'bands': len(arrays2combine),
'band names': arrays2combine,
'data ignore value': nodataVal}
return {'desc': 'masks', 'row_start': 0, 'row_end': self.shape_fullArr[0],
......
......@@ -175,7 +175,14 @@ def L2A_map(L1C_objs, block_size=None, return_tiles=True):
clipextent=clip_extent, clipextent_prj=L2A_obj.arr.prj)
# 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]
if len(L2A_objs) > 1:
L2A_obj = L2A_P.L2A_object.from_sensor_subsystems(L2A_objs)
else:
L2A_obj = L2A_objs[0]
# update attributes
L2A_obj.calc_mask_nodata(overwrite=True)
L2A_obj.calc_corner_positions()
# write output
if CFG.exec_L2AP[1]:
......
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