Extracting the central strip from LANDSAT 7 imagery

Extracting the central strip from LANDSAT 7 imagery

February 8, 2016

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)

 

you may also like:

Presenting Our UAS Research at the UFS Science Team Meeting

Presenting Our UAS Research at the UFS Science Team Meeting

At the recent UFS (Umweltforschungsstation Schneefernerhaus) science team meeting, our group had the opportunity to present current and planned UAS-based research activities in the high-mountain environment surrounding the Zugspitze. The meeting offered an ideal...

BetaFor DFG Project Undergoes Evaluation in Würzburg

BetaFor DFG Project Undergoes Evaluation in Würzburg

In a productive and inspiring atmosphere, our BetaFor DFG project recently welcomed a panel of internationally renowned scientists for its official evaluation. The visit combined the traditional presentation format with a hands-on excursion into the project’s core...

🎄 A Gentle Touch of Christmas in the Foyer

🎄 A Gentle Touch of Christmas in the Foyer

As the year draws to a close and deadlines, reports, and project milestones converge, December often feels like a sprint. Yet this year, something is different the moment you step through the doors of our building. ✨ A Christmas Tree that Tells a Story Standing...

MainPro Project Meeting: Next Project Phase Launched

MainPro Project Meeting: Next Project Phase Launched

On Friday, November 21st, a MainPro project meeting was held to exchange ideas with our cooperating small and medium-sized enterprises (SMEs). Within our MainPro project, we aim to identify the changes caused by climate change in the Main Valley and its surroundings...