Extracting the central strip from LANDSAT 7 imagery

• Tweet

• Tweet

Here is a simple Python code to extract the central strip from Landsat 7 imagery (SLC-off),  that is not affected by the SLC failure. The algorithm shrinks the striping zones through a morphological filter (erosion) and creates a new shapefile AOI that extracts the desired raster extent without striping effects. The code is based on Python for ArcGIS (arcpy) – so you require a ArcGIS license.

General steps:

1. Loop through all Landsat 7 data folders
2. Stack bands for each image
3. Create a mask
4. Erode the mask by 20 pixels
5. Convert the mask to polygon
6. Create a minimum bounding box
7. Clip the original raster through the bbox

import arcpy
from arcpy.sa import *

import sys,os

#  Environment settings (Activate Spatial Analyst, Overwrite Outputs allowed and TIFF compression is LZW)
arcpy.CheckOutExtension("spatial")
arcpy.env.overwriteOutput = True
arcpy.env.compression = 'LZW'

# this is your main directory with all unzipped Landsat datasets
 rootdir = "D:\\DATA\\Landsat7\\"

# create scratch folder "temp" 
temp = "D:\\DATA\\temp\\"

# loop through directory with all unzipped Landsat 7 folders
 for files in os.listdir(rootdir):   
    files = os.path.join(rootdir, files)   
    
    # for each loop the subdir "files" is now the current workspace 
    # (e.g. LE71520322015157-SC20160224113319) that contains the Landsat bands
    arcpy.env.workspace = files  
    rasters = arcpy.ListRasters("*", "TIF")  
    
    # create empty array
    stack_liste = []  
    # loop through all rasters in subdir
    for raster in rasters:   

        image = arcpy.Raster(raster) 
        name  = image.name 
        index = name.split("_")[0]  

        # fill up the array only with the actual spectral bands        
        sr = "_sr_band"  
        if sr in raster:   
            stack_liste.append(raster)             

    # now stack all bands within the array
    stack_name = files + "\\" + index + "_stack.tif"    
    arcpy.CompositeBands_management(stack_liste, stack_name)  

    # convert the image stack to a mask by logical operation with an absurd value that will result in an output "0"
    con = EqualTo(stack_name, 123456789)  

    # now shrink the raster mask with value "0" by 20 pixels
    shrink = temp + "shrink"  
    shrinking = Shrink(con, 20, 0) 
    shrinking.save(shrink)  

    zone = temp + "zone.shp" 
    bbox = temp + "bbox.shp"  

    # conver the shrunk mask to polygon and create a minimum bounding box
    arcpy.RasterToPolygon_conversion(shrink, zone, "NO_SIMPLIFY", "VALUE") 
    arcpy.MinimumBoundingGeometry_management(zone, bbox, "RECTANGLE_BY_WIDTH", "NONE")  

    # now use that bounding box as a mask to cut out the central nadir strip from the original stack
    # Final result 
    extract = files + "\\" + index + "_aoi.tif"  
    ExtractByMask = arcpy.sa.ExtractByMask(stack_name, bbox) 
    ExtractByMask.save(extract)