[GIS] Cloud Removal from Landsat Data

landsatopen-source-gisremote sensing

I'm trying to remove clouds from Landsat 5,7 and 8 data. I'm also trying to do this using only open source options, avoiding ArcGIS. I'm familiar with python and have some QGIS experience. Originally, I was trying to use the Semi-Automatic Classification model in QGIS, but it is giving me a lot of problems, and doesn't always output actual classifications (it often just returns a uniformly valued square…). Can anyone point me in the direction of a better option?

Best Answer

I have managed to identify cloud pixels in Landsat 8 scenes using the following code. I have included some python / pseudo-code below which will return possible cloud pixels. This is the first pass of identifying the potential cloud layer taken from this paper by Zhu and Woodcock (2014).

The paper also details a second pass that further refines the result by excluded misidentified cloud pixels. You should have a look at that too. The method was developed for Landsat 7 so it should work for all Landsat cases. 

def calc_ndsi():
     green = get_green_band()
     swir = get_swir_band()
     return (green - swir)/(green + swir)

def calc_ndvi():
    nir = get_nir_band()
    red = get_red_band()
    return (nir - red)/(nir + red)

def calc_basic_test():
    band_7_test = get_band_seven()
    temp_test = where temp is less than 27 deg C
    ndsi_test = where calc_ndsi < 0.8
    ndvi_test = where calc_ndvi < 0.8
    basic_test = np.logical_and.reduce((band_7_test,
                                    temp_test,
                                    ndsi_test,
                                    ndvi_test))
    return basic_test

def calc_whiteness():
    blue = get blue
    green = get green
    red = get red
    mean_vis = (blue + green + red) / 3
    whiteness = (np.abs((blue - mean_vis)/mean_vis) + 
             np.abs((green - mean_vis)/mean_vis) +
             np.abs((red - mean_vis)/mean_vis))
    whiteness[np.where(whiteness>1)] = 1
    return whiteness

def calc_whiteness_test():
    whiteness_test = where calc_whiteness() < 0.7
    return whiteness_test

def calculate_hot_test():
    band = 'rtoa_'
    blue = get blue
    red = get red
    hot_test = (blue - 0.5*red - 0.08) > 0
    return hot_test

def swirnir_test():
    nir = get NIR
    swir = get SWIR
    return (nir/swir) > 0.75

def calc_pcp():
    return np.logical_and.reduce((calc_basic_test(),
                              calc_whiteness_test(),
                              calculate_hot_test(),
                              swirnir_test()))

So, calc_pcp() will identify the potential cloud pixels in your image.

Note that this method uses the TIRS bands (to get the temperature) in Landsat 8 which are currently not working. Additionally, this method uses data which has already undergone atmospheric correction.

Identifying clouds in remote sensing images is more than a colour processing issue, although calculating the 'whiteness' of the cloud is included in the test.

You should compare the results you get with the Landsat 8 Quality Assessment band.

Related Question