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:
- Loop through all Landsat 7 data folders
- Stack bands for each image
- Create a mask
- Erode the mask by 20 pixels
- Convert the mask to polygon
- Create a minimum bounding box
- 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)






