Módulo 04

Introdução

Usualmente após o download dos arquivos raster, realiza-se recortes das bandas para realizar as análises na área de interesse. Dependendo o número de arquivos que deseja-se fazer o recorte, a atividades pode torna-se cansativa e consumir um tempo significativo. Como trata-se de um processo repetitivo, o objetivo deste post é implementar um script para agilizar o processamento das imagens.

Material e métodos

Entre as possíveis alternativas, uso o método gdal.Warp ( GDAL - Geospatial Data Abstraction Library) para fazer o recorte. A sequência básica para o recorte das imagens, utilizando um arquivo shape, consiste nos seguintes passos:

Crie uma pasta (diretório) e insira na mesma todos os arquivos raster que deseja-se fazer o recorte;
Na mesma pasta, insira o arquivo shape que servirá como máscara, ou seja, fornecerá o contorno do recorte.
Verifique se todos os arquivos estão no mesmo CRS (sistema de referência de coordenadas e.g. EPSG 4326), caso não estejam faça a reprojeção conforme necessário.
Ainda na mesma pasta, crie o arquivo crop_raster.py para inserirmos o código em Python.

O código

Importando as bibliotecas:

from osgeo import ogr
from osgeo import gdal
import glob
import numpy as np

Lendo os arquivos raster contidos na pasta:

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

Se necessário, no código acima, edite a extensão do arquivo raster, e.g. de jp2 para tif. Uma vez criada a lista com todos os arquivos raster que deseja-se recortar, cria-se um objeto com o caminho do arquivo shape que será usado como máscara.

pathShapeMask = str('./shape_mask/shape_mascara.shp')

Antes de realizar o recorte dos arquivos raster, é interessante virificar se os arquivos estão no mesmo CRS. Para verificar o CRS do arquivo shape, uma das opções é:

# Lendo o arquivo shape
ds = ogr.Open(pathShapeMask, 1)  # 0 means read-only. 1 means writeable

# Obtendo informacoes
layer = ds.GetLayer()
spatialRef = layer.GetSpatialRef()

print(spatialRef.ExportToWkt())
>>> 'PROJCS["WGS 84 / UTM zone 21S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-57],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32721"]]'

# fechando o arquivo
ds = None

Para verificar o CRS do arquivo raster:

# Lendo arquivo raster e verificando o CRS
pathRasterFile = './raster/T21LTD_20200305T141041_B02_10m.jp2'
dsR = gdal.Open(pathRasterFile)

print(dsR.GetProjection())
>>> 'PROJCS["WGS 84 / UTM zone 21S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-57],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32721"]]'

dsR = None

Uma vez verificado o CRS, realiza-se o recorte dos arquivos raster. No nome de cada arquivo recortado será adicionado o termo _cut.jp2.

# recortando e renomeando os aquivos raster da lista
for idx in range(len(lsRaster)):
    rasterFile = lsRaster[idx]
    outputFilePath = str('./' + lsRaster[idx].split('.')[1] + '_cut.jp2')
    dsR = gdal.Open(rasterFile)
    dsClip = gdal.Warp(outputFilePath, dsR,
                       cutlineDSName=pathShapeMask,
                       cropToCutline=True, dstNodata=np.nan)
    print("\n File: \n", outputFilePath, " was created \n")
    dsR = None
    dsClip = None

Resultados

Consulte o vídeo no canal do YouTube: ‘YouTube’

Considerações

Os arquivos recortados foram gerados com sucesso. Salienta-se que o desenvolvimento do código pode custar mais tempo que fazer a operação manual para um número reduzido de aquivos. Todavia, uma vez feito o código, o mesmo poderá ser utilizado em outras situações similares.

Outra consideração é que quanto maior o número de aquivos a ser recortado, maior a vantagem do uso do script em relação ao processo manual.

The power of GDAL and OGR

Recorte de arquivos raster com uso de máscara em Python