如何用Python导入不同Landsat平台的指定波段以计算NDVI？

Hey there! I totally get how confusing it can be to figure out which Landsat bands to import when you're just starting out with Python and ArcPy. Let's walk through a solution that'll automatically grab the right red and NIR bands for each of your Landsat 5,7,8 scenes and calculate NDVI for you.

Here's a complete, commented code snippet tailored to your needs:

import os
import arcpy
import re
from arcpy import env
from arcpy.sa import *

# Set your root directory where all unzipped Landsat scene folders are stored
root_dir = r"E:\Thesis\005005-006004\005005"

# Allow overwriting existing files so we don't hit errors during repeated runs
env.overwriteOutput = True

# Check out the Spatial Analyst extension (required for raster calculations)
arcpy.CheckOutExtension("Spatial")

# Loop through each folder in your root directory
for scene_folder in os.listdir(root_dir):
    scene_path = os.path.join(root_dir, scene_folder)
    
    # Skip any items that aren't folders
    if not os.path.isdir(scene_path):
        continue
    
    # Use regex to identify which Landsat mission the scene belongs to
    # Landsat scene IDs start with LT05 (5), LE07 (7), LC08 (8)
    mission_match = re.search(r'L[TE]0(\d)|LC0(\d)', scene_folder)
    if mission_match:
        # Extract the mission number (5,7,8)
        landsat_num = int(mission_match.group(1) or mission_match.group(2))
        
        # Define the correct red and NIR bands based on the mission
        if landsat_num in [5,7]:
            red_band_id = 3
            nir_band_id = 4
        elif landsat_num == 8:
            red_band_id = 4
            nir_band_id = 5
        else:
            print(f"Oops, unsupported Landsat mission: {landsat_num} in folder {scene_folder}")
            continue
        
        # Find the actual TIFF files for our target bands
        red_band_path = None
        nir_band_path = None
        
        for file_name in os.listdir(scene_path):
            # Landsat band files usually end with "_Bx.TIF" (x = band number)
            if file_name.endswith(f"_B{red_band_id}.TIF"):
                red_band_path = os.path.join(scene_path, file_name)
            elif file_name.endswith(f"_B{nir_band_id}.TIF"):
                nir_band_path = os.path.join(scene_path, file_name)
        
        # If we found both bands, proceed to calculate NDVI
        if red_band_path and nir_band_path:
            print(f"Processing scene: {scene_folder}")
            print(f"Red band: {red_band_path}")
            print(f"NIR band: {nir_band_path}")
            
            # Convert the band files to ArcPy Raster objects
            red_raster = Raster(red_band_path)
            nir_raster = Raster(nir_band_path)
            
            # Calculate NDVI using the formula (NIR - Red)/(NIR + Red)
            ndvi_raster = (nir_raster - red_raster) / (nir_raster + red_raster)
            
            # Save the NDVI output to the same scene folder
            ndvi_output_path = os.path.join(scene_path, f"{scene_folder}_NDVI.TIF")
            ndvi_raster.save(ndvi_output_path)
            
            print(f"Success! NDVI saved to: {ndvi_output_path}\n")
        else:
            print(f"Missing red or NIR band in folder {scene_folder}\n")
    else:
        print(f"Couldn't identify Landsat mission for folder {scene_folder}\n")

# Check the extension back in when we're done
arcpy.CheckInExtension("Spatial")

Key Details to Note: #

• Mission Detection: The regex looks for standard Landsat scene ID prefixes (like LT05 for Landsat 5) to automatically determine which bands to use—no manual sorting needed!
• Band Matching: We target files ending with _Bx.TIF (the standard naming convention for Landsat bands) to find the correct red and NIR files.
• Spatial Analyst: The code checks out the Spatial Analyst extension (required for raster math operations like NDVI) and checks it back in when finished.
• Output: NDVI files are saved directly in their respective scene folders with a clear name (e.g., LC08_123456_20200101_NDVI.TIF).

If you run into issues with file naming (e.g., your bands use a different suffix), just adjust the endswith() checks to match your actual file names!

内容的提问来源于stack exchange，提问作者D. Lubenow