Hi,
I am using the BQA layer from Landsat8 to extract a cloud mask via decodeLs8Qa(sample.BQA).cloud.
But somehow the results are very off.
the Landsat image:
The Cloud mask (extracted from the BQA band):
The resulting masked image:
Does anyone have an idea why the clud masking is so off? And how can I improve this?
This is my script:
time_interval = ["2014-01-01", "2016-12-31"]
maxcc = 0.4
EOPATCH_FOLDER = os.path.join(".", "eopatches")
evalscript = """
//VERSION=3
function setup() {
return {
input: [{ bands: ["B01","B02","B03","B04","B05","B06","B07","BQA","dataMask"] }],
output: [
{id:'BANDS', bands: 7, sampleType: 'FLOAT32'},
{id:'cloud_mask', bands: 1, sampleType: 'INT8'},
{id:'cloud_mask_bqa', bands:1, sampleType: 'INT8' },
{id:'mask', bands: 1, sampleType: 'INT8'}
]
}
}
function create_mask(b1,b4,b7){
var abs = Math.abs;
var ceiling = Math.ceil;
var cos = Math.cos;
var exp = Math.exp;
var floor = Math.floor;
var log = Math.log;
var sin = Math.sin;
var sqrt = Math.sqrt;
var truncate = Math.trunc;
var round=Math.round;
var band1 = round(b1*65535);
var band4 = round(b4*65535);
var band7 = round(b7*65535)
let cloud_val;
let output_first = 2.16246741593412 - 0.796409165054949*band4 + 0.971776520302587*sqrt(abs(0.028702220187686*band7*band1 + 0.971297779812314*sin(band1))) + 0.0235599298084993*floor(0.995223926146334*sqrt(abs(0.028702220187686*band7*band1 + 0.971297779812314*sin(band1))) + 0.00477607385366598*abs(0.028702220187686*band7*band1 + 0.971297779812314*sin(band1))) - 0.180030905136552*cos(band4) + 0.0046635498889134*abs(0.028702220187686*band7*band1 + 0.971297779812314*sin(band1));
let output_second = band7;
if (output_first < output_second) {
return(cloud_val = 0);
}
else {
return(cloud_val = 1);
}
}
function evaluatePixel(sample) {
return {
BANDS: [sample.B01, sample.B02, sample.B03, sample.B04, sample.B05, sample.B06, sample.B07, sample.dataMask],
cloud_mask: [create_mask(sample.B01, sample.B04, sample.B07)],
cloud_mask_bqa: [decodeLs8Qa(sample.BQA).cloud],
mask: [sample.dataMask]
}
}
"""
input_task = SentinelHubEvalscriptTask(
features=[(FeatureType.DATA, 'BANDS'),(FeatureType.MASK, 'mask'),(FeatureType.MASK, 'cloud_mask'), (FeatureType.MASK, 'cloud_mask_bqa')],
evalscript=evalscript,
data_collection=DataCollection.LANDSAT_OT_L2, #Landsat OLI-TIRS Level 2
resolution=30,
maxcc=maxcc,
time_difference=datetime.timedelta(hours=2),
config=config
)
band_names = ["B01", "B02", "B03", "B04", "B05", "B06", "B07", "BQA"]
ndvi = NormalizedDifferenceIndexTask(
(FeatureType.DATA, "BANDS"), (FeatureType.DATA, "NDVI"), [band_names.index("B05"), band_names.index("B04")]
)
save = SaveTask(EOPATCH_FOLDER, overwrite_permission=OverwritePermission.OVERWRITE_PATCH)
workflow_nodes = linearly_connect_tasks(input_task, ndvi, save)
workflow = EOWorkflow(workflow_nodes)
for idx in range(2): #TILE_IDS#len(bbox_list)
result = workflow.execute(
{
workflow_nodes[0]: {"bbox": bbox_list[idx], "time_interval": time_interval}, #bbox_list[idx]
workflow_nodes[-1]: {"eopatch_folder": f"landsat8_{idx}"},
}
)
Thanks!
