I'm creating a composite image using Landsat 8 as I want a monthly composite of an area without clouds. I'm having trouble with masking the clouds. I'm using this dataset and I am applying the code below in order to get the images, filter by date, area and cloud cover.
def apply_scale_factors(image):
optical_bands = image.select('SR_B.').multiply(0.0000275).add(-0.2)
thermal_bands = image.select('ST_B.*').multiply(0.00341802).add(149.0)
return image.addBands(optical_bands, None, True).addBands(
thermal_bands, None, True
)
#Bit 1: Dilated Cloud
#Bit 2: Cirrus (high confidence)
#Bit 3: Cloud
#Bit 4: Cloud Shadow
#Bit 5: Snow
def maskClouds(image):
qa = image.select('QA_PIXEL')
cloud_shadow_bit_mask = (1 << 3)
cloud_bit_mask = (1<< 4)
cloud_mask = qa.bitwiseAnd(cloud_shadow_bit_mask).eq(0).And(qa.bitwiseAnd(cloud_bit_mask).eq(0))
return image.updateMask(cloud_mask)
def filterLandsat8(start_date,end_date,aoi,cloud_filter):
return ee.ImageCollection('LANDSAT/LC08/C02/T1_L2').filterDate(ee.Date(start_date),ee.Date(end_date)).filterBounds(aoi).filter(ee.Filter.lte('CLOUD_COVER',cloud_filter)).map(apply_scale_factors).map(maskClouds)
The problem is that the mask is removing a lot of pixels that are not cloud pixels. I even compared the same image collection without and with the cloud mask and as you can see, the maskClouds function removed data that should have not been removed.
Is something in my code wrong? Is this a problem with the QA_PIXEL band?
EDIT:
This is the code I am running:
def apply_scale_factors(image):
optical_bands = image.select('SR_B.').multiply(0.0000275).add(-0.2)
thermal_bands = image.select('ST_B.*').multiply(0.00341802).add(149.0)
return image.addBands(optical_bands, None, True).addBands(
thermal_bands, None, True
)
#Bit 1: Dilated Cloud
#Bit 2: Cirrus (high confidence)
#Bit 3: Cloud
#Bit 4: Cloud Shadow
#Bit 5: Snow
def maskClouds(image):
qa = image.select('QA_PIXEL')
cloud_shadow_bit_mask = (1 << 3)
cloud_bit_mask = (1<< 4)
cloud_mask = qa.bitwiseAnd(cloud_shadow_bit_mask).eq(0).And(qa.bitwiseAnd(cloud_bit_mask).eq(0))
return image.updateMask(cloud_mask)
def filterLandsat8(start_date,end_date,aoi,cloud_filter):
return ee.ImageCollection('LANDSAT/LC08/C02/T1_L2').filterDate(ee.Date(start_date),ee.Date(end_date)).filterBounds(aoi).filter(ee.Filter.lte('CLOUD_COVER',cloud_filter)).map(apply_scale_factors).map(maskClouds)
START_DATE = '2023-10-01'
END_DATE = '2023-10-31'
CLOUD_FILTER = 40
caop = ee.FeatureCollection('projects/ee-rutesantosceiia/assets/caop')
porto_caop = caop.filter(ee.Filter.eq('Distrito','Porto'))
AOI = ee.Geometry(porto_caop.geometry())
#get image collection
landsat8October23Collection = filterLandsat8(START_DATE,END_DATE,AOI,CLOUD_FILTER)
#composition for month of October
landsat8OriginalMedian = landsat8October23Collection.median().clip(AOI)
landsat8OriginalMedianCloudFree = landsat8October23Collection.map(maskClouds).median().clip(AOI)
#run maps
center = AOI.centroid(10).coordinates().reverse().getInfo()
rgb_l8_vis = {
'bands': ['SR_B4', 'SR_B3', 'SR_B2'],
'min': 0.0,
'max': 0.3,
}
labels = [
'Median Landsat 8 Composite for October 2023',
'Median Cloud-Free Landsat 8 Composite for October 2023',
]
geemap.linked_maps(
rows=1,
cols=2,
height="400px",
center=center,
zoom=10,
ee_objects=[landsat8OriginalMedian,landsat8OriginalMedianCloudFree],
vis_params=[rgb_l8_vis,rgb_l8_vis],
labels=labels,
label_position="topright",
)
I have tried different approaches to this but if I mask the cloud bit (3) it always removes parts of the image that I cannot detect any clouds. I will leave a better image of the composite without the cloud mask and with the cloud mask.
Best Answer
The QA_PIXEL in those regions is set at "QA_PIXEL: 55052". The Landsat user's guide says that's "cirrus cloud with high confidence".
So I'm guessing the aerosol band used to do atmospheric correction is saying those areas contain cirrus clouds, which inspecting the RGB data closely, looks to be true: