[GIS] Is it possible to read in specific bands from a multi-band raster with gdal or rasterio

gdalmulti-bandpythonrasterio

I wasn't able to find an answer to this after searching the help files and the forums.

I have a (virtual) raster with 60 bands, where each band represents a realization in time (think of frames of a video). I want to be able to load in cropped, specific band sets from this virtual raster without loading all of them–something like

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

which would create an array with only bands 2,4,6,14, and 35.

I haven't been able to find a way to do this without loading each band separately and filling in a 3D matrix.

Another alternative that I also want to avoid is to re-create each of the individual cropped bands (which are also virtual), then create a new virtual raster stack using only the bands I want, and finally read it into an array.

I saw this was listed as a feature to add to rasterio but I couldn't find any evidence of it actually being added.

Edit: After receiving an answer, I think I should clarify that I'm asking if this can be performed in a single call, rather than implementing a loop to read each band individually.

Best Answer

You can read specific bands in a single call using rasterio by passing a list/tuple of band numbers (Following the GDAL convention, bands are indexed from 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)

Related Solutions

[GIS] Removing some bands or layers from a tiff using gdal

You could create a GDALDataset with as many bands as you have raster bands, then copy the data from each of your bands into the corresponding band in the GDALDataset. Here's some example code in C++ (since that's where I'm most familiar with GDAL).

//create the dataset
const char *filename = "example.tif";
GDALDriver *pDriverTiff = GetGDALDriverManager->getDriverByName("GTiff");
GDALDataset *pRasterDS = pDriverTiff->create(filename, 300, 300, 25, GDT_Float32, NULL);//for create(filename, x size, ysize, number of bands, data type, options)

//fill each raster band
int bands = 25; //this would be the number of bands you want to copy i.e. total bands-bitmaps
int rows = 300; //rows in the raster or 3d array
int cols = 300; //cols in the raster or 3d array

for (int i=0; i<band; i++)
{
    for (int j=0; j<rows; j++)
    {
         for (int k=0; k<cols; k++)
         {
              //here is where you put the code to copy from your 3d array
              // copy array position [i,j,k] to band i, row j, col k, in the GDALDataset
         }

    }   
}

There might be a more efficient way to do this that I am not aware of, but I have used this method many times to add and delete bands from raster datasets.

[GIS] How to create an rgb GeoTIFF file raster from bands using the GDAL python module

You can set the colour interpretation for most rasters using Band.SetColorInterpretation, i.e

outband.SetColorInterpretation(gdal.GCI_RedBand)

However, this is only understood by GDAL based software. To specifically create a GeoTIFF with an RGB baseline PhotometricInterpretation tag (this is not a GDAL tag, but part of the TIFF spec) pass a PHOTOMETRIC=RGB option to the GeoTIFF driver when creating the raster:

driver = gdal.GetDriverByName('GeoTIFF')
options = ['PHOTOMETRIC=RGB', 'PROFILE=GeoTIFF']
outRaster = driver.Create(filepath, cols, rows, bands, datatype, options=options)

You're also trying to write 3 bands to a single band:

(outRaster.GetRasterBand(1).WriteArray(etc...))

Get used to having your arrays in [bands, rows, cols] order not [cols, rows, bands] as this is how gdal arranges them when reading into a numpy array. Then you need to loop through the array "bands" and write each of them to the raster bands:

bands = array.shape[0]
for band in range(bands):
    outRaster.GetRasterBand(band+1).WriteArray( array[band, :, :] )

Try something like:

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

def array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array,format='GTiff'):

    bands = array.shape[0]
    rows = array.shape[1]
    cols = array.shape[2]
    originX = rasterOrigin[0]
    originY = rasterOrigin[1]
    driver = gdal.GetDriverByName(format)
    #Here is where I assign three bands to the raster with int 3
    options = ['PHOTOMETRIC=RGB', 'PROFILE=GeoTIFF']
    outRaster = driver.Create(newRasterfn, cols, rows, bands, gdal.GDT_UInt16, options=options)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    #outband = outRaster.GetRasterBand(1)
    #outband.WriteArray(array)
    for band in range(bands):
        outRaster.GetRasterBand(band+1).WriteArray( array[band, :, :] )

    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(6962)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())

    return outRaster


def example_rgb_creation():
    source = gdal.Open("C:\\RGBtest.tif")
    rgb_array = source.ReadAsArray()

    return rgb_array

exampleRGB = example_rgb_creation()
exampleOutput = "C:\\Temp\\ExampleOutput.tif"
rasterOrigin=[0,0]

outRaster = array2raster(exampleOutput,rasterOrigin,1, 1,exampleRGB)

outRaster.GetRasterBand(1).SetColorInterpretation(gdal.GCI_RedBand)
outRaster.GetRasterBand(2).SetColorInterpretation(gdal.GCI_GreenBand)
outRaster.GetRasterBand(3).SetColorInterpretation(gdal.GCI_BlueBand)

del outRaster

Best Answer

Related Solutions

[GIS] Removing some bands or layers from a tiff using gdal

[GIS] How to create an rgb GeoTIFF file raster from bands using the GDAL python module

Related Question