1

How can I filter areas of scene overlap so that each pixel only has data from one landsat scene? My goal is to have as many images as possible while also maintaining a uniform number of images across the area of interest.

Here is a map showing image count in the collection (darker areas have more images): enter image description here

I have made attempts to filter to only 1 lansdat scene per pixel but can only get something like the following: enter image description here

Here is the code I have been using:


def get_landsat_collection(landsat_num, start_date, end_date, aoi, min_qual, max_cc):
    if landsat_num == 4:
        ls_col = ee.ImageCollection('LANDSAT/LT04/C01/T1_SR').select(
            ['B1', 'B2', 'B3', 'B4', 'B5', 'B7', 'pixel_qa'],
            ['blue', 'green', 'red', 'nir', 'swir1', 'swir2', 'pixel_qa']
        ).filterMetadata(
            'IMAGE_QUALITY', 'not_less_than', min_qual
        ).filterMetadata(
            'CLOUD_COVER', 'not_greater_than', max_cc
        )
    elif landsat_num == 5:
        ls_col = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR').select(
            ['B1', 'B2', 'B3', 'B4', 'B5', 'B7', 'pixel_qa'],
            ['blue', 'green', 'red', 'nir', 'swir1', 'swir2', 'pixel_qa']
        ).filterMetadata(
            'IMAGE_QUALITY', 'not_less_than', min_qual
        ).filterMetadata(
            'CLOUD_COVER', 'not_greater_than', max_cc
        )
    elif landsat_num == 7:
        ls_col = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR').select(
            ['B1', 'B2', 'B3', 'B4', 'B5', 'B7', 'pixel_qa'],
            ['blue', 'green', 'red', 'nir', 'swir1', 'swir2', 'pixel_qa']
        ).filterMetadata(
            'IMAGE_QUALITY', 'not_less_than', min_qual
        ).filterMetadata(
            'CLOUD_COVER', 'not_greater_than', max_cc
        )
    elif landsat_num == 8:
        ls_col = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR').select(
            ['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'pixel_qa'],
            ['blue', 'green', 'red', 'nir', 'swir1', 'swir2', 'pixel_qa']
        ).filterMetadata(
            'IMAGE_QUALITY_OLI', 'not_less_than', min_qual
        ).filterMetadata(
            'CLOUD_COVER', 'not_greater_than', max_cc
        )
    else:
        raise ValueError("Landsat 4, 5, 7, 8 are acceptable options")

    ls_col = ls_col.filterBounds(aoi).filterDate(start_date, end_date)
    return ls_col

def mask_clouds(image):
    qa = image.select('pixel_qa')
    shadow_bit_mask = 1 << 3
    cloud_bit_mask = 1 << 5
    cloud_mask = qa.bitwiseAnd(cloud_bit_mask).eq(0).And(
        qa.bitwiseAnd(shadow_bit_mask).eq(0))
    masked_img = image.updateMask(cloud_mask).copyProperties(image, ['system:time_start'])
    return masked_img

def mask_water(image):
    qa = image.select('pixel_qa')
    water_bit_mask = 1 << 2
    water_mask = qa.bitwiseAnd(water_bit_mask).eq(0)
    masked_img = image.updateMask(water_mask).copyProperties(image, ['system:time_start'])
    return masked_img

def mask_snow(image):
    qa = image.select('pixel_qa')
    snow_mask = 1 << 4
    mask = qa.bitwiseAnd(snow_mask).eq(0)
    masked_img = image.updateMask(mask).copyProperties(image, ["system:time_start"])
    return masked_img

start_yr = 1985
end_yr = 2021

image_qual_min = 9
cloud_cc_max = 50
nys_coords = [
    [-79.99688779, 40.40498417],
    [-71.65014432, 40.40498417],
    [-71.65014432, 44.97328983],
    [-79.99688779, 44.97328983]
]

aoi = ee.Geometry.Polygon(nys_coords)
aoi_ls_collections = []

for year in range(start_yr, end_yr):
    ls_4_col = utils.get_landsat_collection(4, str(year) + '-01-01', str(year) + '-12-31', aoi, image_qual_min,
                                            cloud_cc_max)
    ls_5_col = utils.get_landsat_collection(5, str(year) + '-01-01', str(year) + '-12-31', aoi, image_qual_min,
                                            cloud_cc_max)
    ls_8_col = utils.get_landsat_collection(8, str(year) + '-01-01', str(year) + '-12-31', aoi, image_qual_min,
                                            cloud_cc_max)
    ls_7_col = utils.get_landsat_collection(7, str(year) + '-01-01', str(year) + '-12-31', aoi, image_qual_min,
                                                cloud_cc_max)
    ls_col = ls_4_col.merge(ls_5_col.merge(ls_7_col.merge(ls_8_col)).sort('system:time_start')
    ls_col = ls_col.map(utils.add_ndvi)
    ls_col = ls_col.map(utils.mask_clouds)
    ls_col = ls_col.map(utils.mask_snow)
    ls_col = ls_col.map(utils.mask_water)
    ls_col = ls_col.map(lambda img: img.clip(aoi))
    aoi_ls_collections.append(ls_col)

aoi_ls_collection = aoi_ls_collections.pop(0)
for col in aoi_ls_collections:
    aoi_ls_collection = aoi_ls_collection.merge(col)

input_rp_image = aoi_ls_collection.map(lambda x: ee.Image(ee.Number(x.get("WRS_PATH")).add(ee.Number(x.get("WRS_ROW")).multiply(1000))).toInt16().updateMask(x.select('NDVI').mask()))
input_rp_stack = input_rp_image.toArray().arrayProject([0])
rp_mask = input_rp_stack.arraySort(input_rp_stack).arraySlice(0, -1,).arrayFlatten([["mask"]])
test_collection_masked = aoi_ls_collection.map(lambda x: x.updateMask(ee.Image(ee.Number(x.get("WRS_PATH")).add(ee.Number(x.get("WRS_ROW")).multiply(1000))).toFloat().eq(rp_mask)))

args = dict(
    folder="GEE_DUMP",
    scale=30,
    crs='EPSG:26918',
    region=aoi,
    maxPixels=500000000
)

nys_count_task = geetools.batch.Export.image.toDrive(
    image=test_collection_masked.select('NDVI').count(),
    description="ccdc_filtered_input_count_nys_12142020",
    **args
)

nys_count_task.start()
0

https://developers.google.com/earth-engine/tutorials/tutorial_api_05 read this, should be helping.

Use different reducers to composite images.

3
  • 1
    Thank you for your answer, but that link takes me to a 404 page
    – Lucas
    Dec 14 '20 at 17:27
  • 1
    Kindly, check the URL before sharing :)
    – Yogi
    Dec 14 '20 at 17:33
  • Ok, the link works now. But it doesn't quite solve my problem. The examples there show how to reduce the image collection to one image through mosaicking and compositing. I need each pixel to hold a time series (> 1 observations), and actually the more observations per pixel the better, but I need these observations to be relatively uniform across the area of interest. So for my needs compositing is too aggressive.
    – Lucas
    Dec 15 '20 at 18:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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