Las imágenes satelitales nos han proporcionado la capacidad de ver la superficie de la tierra en los años recientes pero no hemos tenido mucho éxito en entender la dinámica del uso de suelo y su interacción con los factores económicos, sociológicos y políticos. Algunas deficiencias para este este análisis fueron encontradas en el uso de software de GIS comercial pero otras limitaciones están en la forma en que aplicamos procesos lógicos y matemáticos a un set de datos de imágenes satelitales. El manejo de datos geoespaciales en Python nos da un uso eficiente del poder computacional y provee un mayor panorama en análisis de datos.

Este tutorial muestra el procedimiento completo para crear un raster de cambio de uso de suelo proveniente de una comparación entre raster de índice de vegetacion (NDVI) con el uso de Python y las librerías GDAL y Numpy. Se generaron contornos de cambio de uso de suelo con herramientas de GDAL y Osgeo y se realizó un análisis de deforestación basado en los datos de salida y un recuento de imágenes históricas de Google Earth.

Tutorial

Código en Python

Este el el código en Python para este tutorial:

Importar las librerías requeridas

Este tutorial usa solo librerías de Python, Scipy (Numpy, Matplotlib y otros) y GDAL. Para usuarios de Windows, la forma más efectiva para descargar e instalar GDAL es descarga el GDAL Wheel de https://www.lfd.uci.edu/~gohlke/pythonlibs/ e instalarlo desde Anaconda Prompt.

from osgeo import ogr, gdal, osr
import numpy as np
import os
import matplotlib.pyplot as plt

Configurando las rutas de los archivos de entrada y salida

Declaramos las rutas para las bandas de las imágenes Landsat 8 y para los raster de salida: NDVI y cambio de uso de suelo. También se define una ruta para los contornos de salida.

#Input Raster and Vector Paths

#Image-2019
path_B5_2019="../Image20190203clip/LC08_L1TP_029047_20190203_20190206_01_T1_B5_clip.TIF"
path_B4_2019="../Image20190203clip/LC08_L1TP_029047_20190203_20190206_01_T1_B4_clip.TIF"

#Image-2014
path_B5_2014="../Image20140205clip/LC08_L1TP_029047_20140205_20170307_01_T1_B5_clip.TIF"
path_B4_2014="../Image20140205clip/LC08_L1TP_029047_20140205_20170307_01_T1_B4_clip.TIF"

#Output Files

#Output NDVI Rasters 
path_NDVI_2019 = '../Output/NDVI2019.tif'
path_NDVI_2014 = '../Output/NDVI2014.tif'
path_NDVIChange_19_14 = '../Output/NDVIChange_19_14.tif'

#NDVI Contours
contours_NDVIChange_19_14 = '../Output/NDVIChange_19_14.shp'

Abrir las bandas de Landsat 8 con GDAL

En esta parte la banda de rojo (Banda 4) y la del infrarojo cercano (Banda 5) son leídas con comandos de la librería de GDAL y el arreglo de pixeles es leído como matriz con elementos de numero real de 32 bits.

#Open raster bands
B5_2019 = gdal.Open(path_B5_2019)
B4_2019 = gdal.Open(path_B4_2019)
B5_2014 = gdal.Open(path_B5_2014)
B4_2014 = gdal.Open(path_B4_2014)

#Read bands as matrix arrays
B52019_Data = B5_2019.GetRasterBand(1).ReadAsArray().astype(np.float32)
B42019_Data = B4_2019.GetRasterBand(1).ReadAsArray().astype(np.float32)
B52014_Data = B5_2014.GetRasterBand(1).ReadAsArray().astype(np.float32)
B42014_Data = B4_2014.GetRasterBand(1).ReadAsArray().astype(np.float32)

Comparación de los tamaños de matrix y parámetros de geotransformación

La operación entre imágenes de Landsat 8 depende que la resolución, tamaño, sistema de coordenadas y parametros de geotransformación sean los mismos para ambas imágenes. En caso que estas características con coincidan, se necesitará un proceso geoespacial como combar, cortar, reproyectar, escalar, etc.

print(B5_2014.GetProjection()[:80])
print(B5_2019.GetProjection()[:80])
if B5_2014.GetProjection()[:80]==B5_2019.GetProjection()[:80]: print('SRC OK')

PROJCS["WGS 84 / UTM zone 13N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84
PROJCS["WGS 84 / UTM zone 13N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84
SRC OK

print(B52014_Data.shape)
print(B52019_Data.shape)
if B52014_Data.shape==B52019_Data.shape: print('Array Size OK')

(610, 597)
(610, 597)
Array Size OK

print(B5_2014.GetGeoTransform())
print(B5_2019.GetGeoTransform())
if B5_2014.GetGeoTransform()==B5_2019.GetGeoTransform(): print('Geotransformation OK')

(652500.0, 29.98324958123953, 0.0, 2166000.0, 0.0, -30.0)
(652500.0, 29.98324958123953, 0.0, 2166000.0, 0.0, -30.0)
Geotransformation OK

Obtener parámetros de geotransformación

Debido a que las bandas raster de Landsat 8 tienen el mismo tamaño, están en la misma posición y poseen el mismo sistema de coordenadas se pueden obtener parámetros espaciales comunes.

geotransform = B5_2014.GetGeoTransform()

originX,pixelWidth,empty,finalY,empty2,pixelHeight=geotransform
cols =  B5_2014.RasterXSize
rows =  B5_2014.RasterYSize

projection = B5_2014.GetProjection()

finalX = originX + pixelWidth * cols
originY = finalY + pixelHeight * rows

Calcular el NDVI y almacenarlo en un archivo aparte

Podemos calcular el NDVI con las herramientas de álgebra de matrices de Numpy. El resultado puede ser exportado como un raster usando GDAL con los parámetros de geotransformación obtenidos. Para hacer más efectivo y claro el trabajo se ha generado una función para exportar los rasters como archivos Gtiff.

ndvi2014 = np.divide(B52014_Data - B42014_Data, B52014_Data+ B42014_Data,where=(B52014_Data - B42014_Data)!=0)
ndvi2014[ndvi2014 == 0] = -999

ndvi2019 = np.divide(B52019_Data - B42019_Data, B52019_Data+ B42019_Data,where=(B52019_Data - B42019_Data)!=0)
ndvi2019[ndvi2019 == 0] = -999

def saveRaster(dataset,datasetPath,cols,rows,projection):
    rasterSet = gdal.GetDriverByName('GTiff').Create(datasetPath, cols, rows,1,gdal.GDT_Float32)
    rasterSet.SetProjection(projection)
    rasterSet.SetGeoTransform(geotransform)
    rasterSet.GetRasterBand(1).WriteArray(dataset)
    rasterSet.GetRasterBand(1).SetNoDataValue(-999)
    rasterSet = None

saveRaster(ndvi2014,path_NDVI_2014,cols,rows,projection)

saveRaster(ndvi2019,path_NDVI_2019,cols,rows,projection)

Mostrar las imágenes NDVI

Creamos una función para representar los NDVI resultantes junto a una escala de color.

extentArray = [originX,finalX,originY,finalY]
def plotNDVI(ndviImage,extentArray,vmin,cmap):
    ndvi = gdal.Open(ndviImage)
    ds2019 = ndvi.ReadAsArray()
    plt.figure(figsize=(20,15))
    im = plt.imshow(ds2019, vmin=vmin, cmap=cmap, extent=extentArray)#
    plt.colorbar(im, fraction=0.015)
    plt.xlabel('Este')
    plt.ylabel('Norte')
    plt.show()

plotNDVI(path_NDVI_2014,extentArray,0,'YlGn')

plotNDVI(path_NDVI_2019,extentArray,0,'YlGn')

Creamos una imagen de cambio de uso de suelo

Se crea una imagen de cambio de uso de suelo basado en la diferencia de los NDVI del 2019 con el 2014. La imagen es exportada como un archivo Gtiff.

ndviChange = ndvi2019-ndvi2014
ndviChange = np.where((ndvi2014>-999) & (ndvi2019>-999),ndviChange,-999)
ndviChange

array([[-9.9900000e+02,  3.5470784e-02,  3.7951261e-02, ...,
         3.1590372e-02,  3.8002759e-02,  2.6692629e-02],
       [-9.9900000e+02,  3.7654012e-02,  5.8923483e-02, ...,
        -5.8691800e-03,  1.8922418e-02,  1.9506305e-02],
       [-9.9900000e+02,  3.2184020e-02,  3.7395060e-02, ...,
        -6.3773081e-02, -3.0675709e-02,  3.8942695e-02],
       ...,
       [-9.9900000e+02, -1.6597062e-02,  1.1402190e-02, ...,
         2.7218461e-03,  2.5526762e-02,  6.6639006e-02],
       [-9.9900000e+02,  5.9944689e-03, -4.0770471e-03, ...,
         5.7501763e-02,  4.6757758e-03,  8.4484324e-02],
       [-9.9900000e+02, -9.9900000e+02, -9.9900000e+02, ...,
        -9.9900000e+02, -9.9900000e+02, -9.9900000e+02]], dtype=float32)

saveRaster(ndviChange,path_NDVIChange_19_14,cols,rows,projection)

plotNDVI(path_NDVIChange_19_14,extentArray,-0.5,'Spectral')

Creamos contornos

Finalmente se generan las líneas de contorno de los valores de NDVI. Esto se hace también con herramientas de la librerí GDAL y Osgeo.

Dataset_ndvi = gdal.Open(path_NDVIChange_19_14)#path_NDVI_2014
ndvi_raster = Dataset_ndvi.GetRasterBand(1)

ogr_ds = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(contours_NDVIChange_19_14)

prj=Dataset_ndvi.GetProjectionRef()#GetProjection()

srs = osr.SpatialReference(wkt=prj)#
#srs.ImportFromProj4(prj)

contour_shp = ogr_ds.CreateLayer('contour', srs)
field_defn = ogr.FieldDefn("ID", ogr.OFTInteger)
contour_shp.CreateField(field_defn)
field_defn = ogr.FieldDefn("ndviChange", ogr.OFTReal)
contour_shp.CreateField(field_defn)
#Generate Contourlines
gdal.ContourGenerate(ndvi_raster, 0.1, 0, [], 1, -999, contour_shp, 0, 1)
ogr_ds = None