0

I followed the instructions here. My input layers are 7 .tif images, each representing a spectral band of a Landsat image. I'm attempting to extract these bands using GDAL:

raster_dataset = gdal.Open(raster_data_path, gdal.GA_ReadOnly)

bands_data = []
for b in range(1, raster_dataset.RasterCount+1):
    band = raster_dataset.GetRasterBand(b)
    bands_data.append(band.ReadAsArray())

Issue is, raster_dataset.RasterCount == 1 and I'm getting arrays full of zeros.

mockitodorito
  • 1
  • 1
  • 1
  • If you have 7 TIFF images for one scene each will have 1 band, therefore 7 raster datasets to read. Do you have enough memory to read the entire image as an array? It might be better to use ReadRaster to read a block than ReadAsArray to try to grab the entire band as one array. – Michael Stimson Mar 21 '17 at 03:36
  • @MichaelMiles-Stimson Not sure about the memory question - haven't gotten that far. I'm trying to follow this guy's blog post: https://www.machinalis.com/blog/python-for-geospatial-data-processing/#training. He seems to do it just fine and then transforms the band’s data into a numpy array. – mockitodorito Mar 21 '17 at 03:47
  • Regardless, I'm wondering if the image was somehow merged incorrectly since raster_dataset.RasterCount == 1 and not 7. – mockitodorito Mar 21 '17 at 03:50
  • Reading a whole image less than a megabyte works well but doesn't work so well for huge scenes, especially if you exceed your addressable memory. What are you trying to achieve here? Is the problem that the image doesn't appear to be stacked? To create a stacked image you need to create it with the correct number of bands, did you do that in python or in QGIS toolbox?. – Michael Stimson Mar 21 '17 at 03:51
  • @MichaelMiles-Stimson Correct. The image doesn't appear to be stacked. Shouldn't the instructions I linked in my original post do that? I'm not sure what step I missed. – mockitodorito Mar 21 '17 at 03:52
  • I'm not sure, I've not used that tool, I can't see where the band count is modified but that should be calculated from the input bands. It appears that the image stacking isn't successful if the band count is still 1. Is the data type correct? the image shows Float32 but your landsat should be Unsigned Char or UInt16. Were there any error messages when you ran the stacking tool? – Michael Stimson Mar 21 '17 at 03:59
  • Perhaps something like this http://gis.stackexchange.com/questions/193958/image-stacking-with-python-gdal might work better for you. No need to convert to a single TIFF if you don't need to, a VRT will allow you to read the native images as a band of a composite. Or http://gis.stackexchange.com/questions/22300/layer-stacking-rasters-in-qgis – Michael Stimson Mar 21 '17 at 04:07

1 Answers1

2

I do it exactly as you've written it, here's my production code:

for band_n in range(n_parameters):
    band = ds.GetRasterBand(band_n + 1)
    band_array = band.ReadAsArray()
    all_data[band_n] = band_array

So I don't think there's an issue with this code.

I'd check your source files in something like QGIS, where you can view the data and ensure that it's correct in the files.

Also, you mention you have 7 files, so don't you need to loop over each file rather than through bands in one file?

Alex Leith
  • 13,453
  • 30
  • 69