Como cortar múltiples bandas de una imagen Landsat 8 con Python y GDAL - Tutorial

Los estudios del ciclo del agua, uso de suelo y ecosistemas requieren de la evaluación de imágenes satelitales de largos periodos de observación. Las imágenes satelitales necesitan ser recortadas a un área de interés, filtradas, escaladas y convertidas a un formato deseado y almacenado en un archivo. Actualmente, las aplicaciones de GIS son totalmente capaces para el manejo y análisis espacial cuando la cantidad de imágenes ráster es limitada; sin embargo, cuando se trabaja con grandes grupos de imágenes el procesamiento espacial en una interfaz gráfica puede ser lenta y poco práctico. El uso de lenguajes de programación y procesamiento como Python y avanzadas librerias espaciales como GDAL (gdal.org) ayudan a la transformación de datos espaciales de un modo más abstracto y eficaz. Este tutorial muestra el procedimiento completo para cortar un set de bandas de una imagen Landsat 8 y almacenarlo en un archivo renombrado con un sufijo de cada banda.

Este tutorial se realiza en una plataforma de programación para Python llamado Jupyter Notebook. Los archivos de entrada que son bandas raster y shapefile del área de interés (AOI) tienen que estar en un mismo sistema de referencia de coordenadas (SRC) de lo contrario la librería GDAL no localizará los datos espaciales en la posición correcta. Este tutorial muestra el procedimiento para los sets de bandas de una imagen Landsat 8, un ejemplo para una única banda que se proporciona en los scripts de los datos de entrada. Finalmente, este tutorial muestra el raster completo y cortado en un software de Sistemas Información Geográfica como QGIS.

Para los usuarios principiantes GIS con poca expedicionaria de programación, se aconseja terminar el tutorial con los datos proporcionados en la parte de datos de entrada de este tutorial. Una vez que el usuario tenga más experiencia con el software, el usuario puede modificar el código para el procesamiento personalizado de imagenes ráster.

Tutorial

Para la instalación de GDAL y otras librerías espaciales por favor revisar este tutorial.

Código Python

Este es el código Python para cortar bandas múltiples de una imágen Landsat 8:

Importar las librerías requeridas

Este tutorial solo requiere de una librería GDAL y las librerías de Python core como la librería Operating System (OS).

from osgeo import gdal
import os

Importar las librerías requeridas

Este tutorial solo requiere de una librería GDAL y las librerías de Python core como la librería Operating System (OS).

#Input and output paths
#inputPath = '../Images/'
inputPath = '../Images_c1/'
outputPath = '../Output/'

#Input Raster and Vector Paths
bandList = [band for band in os.listdir(inputPath) if band[-4:]=='.TIF']
bandList

['LC08_L1TP_029047_20190203_20190206_01_T1_B10_c1.TIF',
 'LC08_L1TP_029047_20190203_20190206_01_T1_B11_c1.TIF',
 'LC08_L1TP_029047_20190203_20190206_01_T1_B1_c1.TIF',
 'LC08_L1TP_029047_20190203_20190206_01_T1_B2_c1.TIF',
 'LC08_L1TP_029047_20190203_20190206_01_T1_B3_c1.TIF',
 'LC08_L1TP_029047_20190203_20190206_01_T1_B4_c1.TIF',
 'LC08_L1TP_029047_20190203_20190206_01_T1_B5_c1.TIF',
 'LC08_L1TP_029047_20190203_20190206_01_T1_B6_c1.TIF',
 'LC08_L1TP_029047_20190203_20190206_01_T1_B7_c1.TIF',
 'LC08_L1TP_029047_20190203_20190206_01_T1_B8_c1.TIF',
 'LC08_L1TP_029047_20190203_20190206_01_T1_B9_c1.TIF',
 'LC08_L1TP_029047_20190203_20190206_01_T1_BQA_c1.TIF']

#Shapefile of Area of Influence
shp_clip = '../Shp/AOI_1.shp'

Modificación de raster

Cortar el raster seleccionado con la opción Warp de GDAL.

for band in bandList:
    print(outputPath + band[:-4]+'_c2'+band[-4:])
    options = gdal.WarpOptions(cutlineDSName=shp_clip,cropToCutline=True)
    outBand = gdal.Warp(srcDSOrSrcDSTab=inputPath + band,
                        destNameOrDestDS=outputPath + band[:-4]+'_c2'+band[-4:],
                        options=options)
    outBand= None

Datos de entrada

Puedes descargar los datos de entrada y scripts del tutorial de este enlace.