The power of GDAL and OGR

Cálculo de NDVI com GDAL e Python

ogr / gdal / python / geoprocessing / NDVI

Módulo 05

Introdução

Um dos índices de vegetação utilizados com frequência é o NDVI (índice de vegetação por diferença normalizada) (KRIEGLER et al., 1969), o qual é oriundo das banda infra-vermelho (NIR) próximo e vermelho (RED):

$$ NDVI = \frac{NIR - RED}{NIR + RED} $$

Há diversos recursos para a realização dos cálculos do NDVI, entre eles cita-se o QGIS, GRASS, códigos próprios em diversas linguagens de programação, etc. Como mais uma alternativa disponível, e para mostrar a utilização de um script (o que deixa o processo mais prática), este post teve o objetivo de demonstrar um script desenvolvido em Python para o cálculo do NDVI por meio da biblioteca GDAL.

Na sequência é apresentado:

  • leitura dos arquivos raster;
  • obtenção de informações dos arquivos raster;
  • Processamento das bandas;
  • Visualização do NFV;
  • Gravação do arquivo em disco.

Como há um vídeo no YouTube (LINK DO VIDEO) explicando o procedimento, não será detalhado o material e métodos e os resultados. Apenas será disponibilizado código:

"""The power of GDAL and OGR: calcular ndvi."""

# ======================================================================
#                                        Rafael Tieppo
#                                        rafaeltieppo@yahoo.com.br
#                                        https://rafatieppo.github.io/
#                                        2022-06-22
# The power of GDAL and OGR: calcular ndvi
# ======================================================================

from osgeo import gdal
import glob
import numpy as np
import matplotlib.pyplot as plt

# ------------------------------------------------------------
# Lendo os arquivos raster contidos na pasta:

lsRaster = glob.glob('./raster/*.jp2')
lsRaster


def find_band(list_lyrs, bandname):
    """Find a layer by name in a list with layers."""
    for idx in list_lyrs:
        if bandname in idx:
            return idx


# ------------------------------------------------------------
# Lendo arquivo raster - b8 iv sentinel

dsB08 = gdal.Open(find_band(lsRaster, 'B08'))
transform = dsB08.GetGeoTransform()
rows = dsB08.RasterYSize
cols = dsB08.RasterXSize
bands = dsB08.RasterCount
prj = dsB08.GetProjection()

print('\n transform: {} \n rows: {} \n cols: {} \n \
 bands: {} \n proj: {} \n'.format(
    transform, rows, cols, bands, prj))

datatype = gdal.GetDataTypeName(dsB08.GetRasterBand(1).DataType)
datatype

arrB08 = dsB08.GetRasterBand(1).ReadAsArray().astype(np.float32)
arrB08

dsB08 = None

plt.imshow(arrB08, cmap='viridis')
plt.colorbar()
plt.show()

# ------------------------------------------------------------
# Lendo arquivo raster b04 red; sentinel
dsB04 = gdal.Open(find_band(lsRaster, 'B04'))
arrB04 = dsB04.GetRasterBand(1).ReadAsArray().astype(np.float32)
dfB04 = None

# ------------------------------------------------------------
# calculo ndvi
ndvi = (arrB08 - arrB04) / (arrB08 + arrB04)
ndvi = np.nan_to_num(ndvi)
ndviRound = np.round(ndvi, 2)

plt.imshow(ndviRound, cmap='viridis')
plt.colorbar()
plt.show()

# ------------------------------------------------------------
# criando e gravando o ndvi como arquivo GTiff

# driver
driver = gdal.GetDriverByName('GTiff')
dsout = driver.Create('./raster/ndvi.tif',
                      xsize=cols, ysize=rows, bands=1,
                      eType=gdal.GDT_Float32,
                      options=["COMPRESS=LZW"])

# definicao de limites e tamanho de pixel
dsout.SetGeoTransform(transform)

# definicao de crs
dsout.SetProjection(prj)

# definindo array
dsout.GetRasterBand(1).WriteArray(ndviRound.astype(np.float32))
# fechando arquivo
dsout = None

Referências

 Share!

 
comments powered by Disqus