Imitating MapBox' Cloudless Atlas

Inspired by several posts which explain MapBox' marvellous Cloudless Atlas, I wanted to imitate a cloudless atlas by myself. Limited in time and computer power I had to do it much simpler of course. But using first class spatial Python libraries I could achieve quite nice results with a simple Python script.

The result: a cloudless Landsat 8 scene

As input images I use all available Landsat 8 images for one scene. I store the image names in an array and define the upper left and lower right corner:

filenames = [
    "LC81280472013103LGN01",
    "LC81280472013119LGN01",
    "LC81280472013151LGN00",
    "LC81280472013167LGN00",
    "LC81280472013279LGN00",
    "LC81280472013295LGN00",
    "LC81280472013327LGN00",
    "LC81280472013343LGN00",
    "LC81280472014010LGN00",
    "LC81280472014026LGN00",
    "LC81280472014058LGN00",
    "LC81280472014090LGN00",
    "LC81280472014106LGN00",
    "LC81280472014154LGN00"
]
   
# upper left corner for LC8128047:
ul = (208800, 2188289)
# lower right corner for LC8128047:
lr = (430215, 1967195)

A small extent in chronological order

Then all input images are read and stored aa numpy array to a list. This is done for each color band:

for file in filenames:
    ds = gdal.Open("%s/%s/%s.jpg" % (datadir, file, file), gdalconst.GA_ReadOnly)

    geotransform = ds.GetGeoTransform()
    projection = ds.GetProjection()
    originX = geotransform[0]
    originY = geotransform[3]
    pixelWidth = geotransform[1]
    pixelHeight = geotransform[5]
        
    ulx = int((ul[0] - originX) / pixelWidth)
    uly = int((ul[1] - originY) / pixelHeight)
    lrx = int((lr[0] - originX) / pixelWidth)
    lry = int((lr[1] - originY) / pixelHeight)
        
    redband = ds.GetRasterBand(1)
    redpixels.append(redband.ReadAsArray(ulx, uly, (lrx - ulx), (lry - uly)))
    greenband = ds.GetRasterBand(2)
    greenpixels.append(greenband.ReadAsArray(ulx, uly, (lrx - ulx), (lry - uly)))
    blueband = ds.GetRasterBand(3)
    bluepixels.append(blueband.ReadAsArray(ulx, uly, (lrx - ulx), (lry - uly)))

I also prepare the output image with the same dimensions as the input images:

height = len(redpixels[0])
width = len(redpixels[0][0])
   
out_redpixels = np.empty((height, width), dtype=int)
out_greenpixels = np.empty((height, width), dtype=int)
out_bluepixels = np.empty((height, width), dtype=int)

Now comes the main part where the pixels from each band and image are sorted. Then the second darkest pixel is selected for the output image. Since the darkest pixels often are cloud shadows, the second or even third darkest pixel can be chosen, depending on the number of available images and cloud cover of the processed scene.

for col in range(width):
    for row in range(height):

        # Choose the second darkest pixel
        j = 1

        out_redpixels[row][col] = np.sort(np.array([rp[row][col] for rp in redpixels]))[j]
        out_greenpixels[row][col] = np.sort(np.array([gp[row][col] for gp in greenpixels]))[j]
        out_bluepixels[row][col] = np.sort(np.array([bp[row][col] for bp in bluepixels]))[j]

After that the output bands are written and georeferenced to the final image using GDAL.

A small extent with pixels ordered by value

In another test I also tried to exclude cloud pixels using the cloud mask and calculating the mean or the median from the remaining pixels, but the result was not satisfying with this limited number of images.