8

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.

Improve this question

2 Answers 2

Reset to default
6

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)
Improve this answer
4

This appears to work with 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)
Improve this answer
4
  • Aye, and you can do the same loop with my original code. I guess I was asking whether it could be done in a single call, rather than reading each band one at a time. I suppose that any solution would be doing pretty much exactly what you've posted, so there wouldn't be a major speedup.
    – Jon
    Jul 5, 2017 at 23:28
  • 1
    That's correct. I'm quite sure this can't be done without a loop or list comprehension, and if it could I doubt there would be any speedup.
    – Loïc Dutrieux
    Jul 5, 2017 at 23:34
  • 3
    Yes you can read multiple bands in a single call, i.e. src.read(bandlist).
    – user2856
    Jul 6, 2017 at 0:25
  • Thanks @Luke, that's very slick!!! rasterio never stops to impress me :)
    – Loïc Dutrieux
    Jul 6, 2017 at 0:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.