3

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. Left Image: Without Cloud Mask,Right Image: With Cloud Mask

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.

without cloud mask with cloud mask

2
  • Please provide a complete example; this code can't be run. What are the dates, roi and cloudy_percentage you used? Commented Nov 15, 2023 at 13:10
  • @NoelGorelick I have added the code I am running. The AOI is the district of Porto in case you cannot get access to the geometry.
    – rrwork___
    Commented Nov 15, 2023 at 14:27

1 Answer 1

1

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:

enter image description here

3
  • Thank you so much. That makes more sense.
    – rrwork___
    Commented Nov 16, 2023 at 16:29
  • Is there an easy way to remove these cirrus clouds without losing the data? I know this could be a completely different question but I was curious if instead of not using this image at all if I could improve it
    – rrwork___
    Commented Nov 16, 2023 at 16:45
  • 1
    Depends on what you're then going to do with the data. (It might be doable for some use cases). See: mdpi.com/2072-4292/9/8/834 Commented Nov 16, 2023 at 20:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.