Es posible leer en bandas específicas de un multi-banda de ráster con gdal o rasterio?

Question

Yo no era capaz de encontrar una respuesta a esta después de buscar en los archivos de ayuda y los foros.

Tengo una (virtual) ráster con 60 bandas, donde cada banda representa una realización en el tiempo (creo que de fotogramas de video). Quiero ser capaz de cargar en recortada, banda específica de los conjuntos de este virtual trama sin la carga de todos ellos, algo como

myVRT = gdal.Open(vrtfile)
bands = myVRT.GetRasterBand([2,4,6,14,35])
cropped_bandstack = bands.ReadAsArray(
   xoff = xoffset,
   yoff = yoffset,
   xsize = xsz,
   ysize = ysz)

que iba a crear una matriz con sólo bandas 2,4,6,14, y 35.

No he sido capaz de encontrar una manera de hacer esto sin la carga de cada banda por separado y el llenado de una matriz 3D.

Otra alternativa que también quiero evitar es volver a crear cada uno de los recortada bandas (que son también virtuales), a continuación, crear un nuevo virtual trama de la pila utilizando sólo las bandas que yo quiero, y finalmente leer en una matriz.

Vi este fue catalogado como una función para añadir a rasterio, pero no pude encontrar ninguna evidencia de que en realidad se añadió.

Edit: Después de recibir una respuesta, creo que debo aclarar que yo estoy preguntando si esto se puede realizar en una sola llamada, en lugar de implementar un bucle para leer cada banda de forma individual.

Answer 1

Usted puede leer específicos de las bandas en una sola llamada utilizando rasterio por el paso de una lista o tupla de la banda de números (Tras la GDAL convención, las bandas están indexadas a partir de 1):

import rasterio
rasterio.__version__
'1.0a8'

dataset = rasterio.open('multiband.tif')

dataset.count
4

dataset.read((1,2)) #read 1st two bands into an array.

array([[[ 85,  98,  75, ...,  53,  55,  55],
        [ 84,  94,  76, ...,  54,  55,  54],
        [ 68,  60,  55, ...,  53,  54,  53],
        ...,
        [ 67,  67,  66, ...,  63,  63,  62],
        [255, 255, 255, ..., 255, 255, 255],
        [255, 255, 255, ..., 255, 255, 255]],

       [[ 78,  88,  65, ...,  41,  43,  45],
        [ 77,  84,  66, ...,  42,  43,  44],
        [ 61,  51,  46, ...,  41,  42,  43],
        ...,
        [ 77,  77,  77, ...,  71,  70,  69],
        [255, 255, 255, ..., 255, 207, 255],
        [255, 255, 255, ..., 191,   0, 135]]], dtype=uint8)

Answer 2

Esto parece que funciona con rasterio

import rasterio
import numpy as np

vstack = '/path/to/virtual_stack.vrt'
bands = [2,4,6,14,35]

# define cropping window ((row_start, row_stop), (col_start, col_stop))
window = ((10, 50), (30, 40))

with rasterio.open(vstack) as src:
    # Create zero array (you may want to set dtype too)
    array = np.zeros((window[0][1] - window[0][0],
                      window[1][1] - window[1][0],
                      len(bands)))
    # Fill the array
    for i, band in enumerate(bands):
        array[:,:,i] = src.read(band, window=window)