Extraer bandas individuales de un ráster RGB con un script de Python

El apilamiento de bandas ráster es fácil de realizar en QGIS con la opción Raster -> Miscellaneous -> Merge -> Layer stack de la barra de menús. Sin embargo, el proceso contrario, es decir, recuperar las bandas individuales a partir del ráster apilado no parece tan obvio. No obstante, los que han trabajado con el plugin de GRASS, probablemente conozcan que la incorporación de un ráster RGB a los datasets de GRASS automáticamente los separa en sus tres bandas originales (no acepta ráster apilados). Esto pudiera ser aprovechado para realizar la individualización de las bandas en un ráster RGB.

Por otra parte, una vía más expedita para la separación de las bandas puede ser realizada a través de un script de Python utilizando la librería GDAL. Para ello, se lee el dataset completo de las bandas apiladas y, para cada banda particular, se crea un ráster individual con su dataset respectivo. Uno de tales procesos pasa por leer los values por bloques completos, almacenarlos como lista, “aplanar” la lista, convertirla en array y luego escribirla como ráster. El código es el siguiente:

from osgeo import gdal, osr
import os, struct
import numpy as np

layer = iface.activeLayer()

provider = layer.dataProvider()

path = provider.dataSourceUri()

fmttypes = {'Byte':'B', 'UInt16':'H', 'Int16':'h', 'UInt32':'I', 'Int32':'i', 'Float32':'f', 'Float64':'d'}

dataset = gdal.Open(path)

#Get projection
prj = dataset.GetProjection()

#setting number band (blue in this case)
number_band = 3

band = dataset.GetRasterBand(number_band)

geotransform = dataset.GetGeoTransform()

# Set name of output raster
if number_band == 3:
    output_file = "/home/zeito/pyqgis_data/blueband.tif"

# Create gtif file with rows and columns from parent raster 
driver = gdal.GetDriverByName("GTiff")

columns, rows = (band.XSize, band.YSize)

print "rows = %d columns = %d" % (columns, rows)

BandType = gdal.GetDataTypeName(band.DataType)

print "Band Type = ", BandType

raster = []

for y in range(band.YSize):

    scanline = band.ReadRaster(0, y, band.XSize, 1, band.XSize, 1, band.DataType)
    values = struct.unpack(fmttypes[BandType] * band.XSize, scanline)
    raster.append(values)

dst_ds = driver.Create(output_file, 
                       columns, 
                       rows, 
                       1, 
                       band.DataType)

#flattened list of raster values
raster = [ item for element in raster for item in element ]

#transforming list in array
raster = np.asarray(np.reshape(raster, (rows,columns)))

##writting output raster
dst_ds.GetRasterBand(1).WriteArray( raster )

#setting extension of output raster
# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
dst_ds.SetGeoTransform(geotransform)

# setting spatial reference of output raster 
srs = osr.SpatialReference(wkt = prj)
dst_ds.SetProjection( srs.ExportToWkt() )

#Close output raster dataset 
dst_ds = None

#Close main raster dataset
dataset = None

Para probar el código se uso el ráster RGB siguiente:

rgb1

La ejecución en la Python Console de QGIS permite obtener:

rgb2

que corresponde a la banda azul del ráster RGB.

Es de hacer notar que este código puede hacerse más sencillo y directo usando el método ‘ReadAsArray’ de GDAL.

Esta entrada fue publicada en GDAL, PyQGIS, QGIS, SIG, Software Libre. Guarda el enlace permanente.

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s