3

I'm trying to read in several rasters as array using gdal in a python script, but I keep getting NoneType returned and can't figure out why.

The process is this:

  1. Use gdalbuildvrt (http://www.gdal.org/gdalbuildvrt.html) to read in several rasters into gdal's virtual format.
  2. Loop through the rasters (aka the layers) of the .vrt, reading them into a temporary array and adding the temp array to the cumulative (or total) array.
  3. Write the total array to the raster.

import os
import gdal
from osgeo.gdal import *
import numpy as np
gdal.UseExceptions() 

# intxtlist = a text file with the rasters to be put in the vrt; in this case there are 3 files
VRTstack = '/path/tempVRT.vrt'
os.system('gdalbuildvrt -separate -vrtnodata 0 -input_file_list %s %s' % (intxtlist, VRTstack))

vrtDS = gdal.Open(VRTstack)
Nlayers = vrtDS.RasterCount # Nlayers = 3 because there were 3 input rasters
Nsamp = vrtDS.RasterXSize
Nline = vrtDS.RasterYSize

totalArr = np.zeros((Nline, Nsamp), dtype = 'uint16') # preallocate an array for the final output total
for N in range(0, Nlayers):
    tempArr = vrtDS.GetRasterBand(N+1).ReadAsArray() # tempArray is now Nth tiff of stack as array; N+1 b/c gdal starts at 1
    print type(tempArr)
    totalArr += tempArr  # add the temp array to the final array
## code here to write totalArr as raster

The problem I'm getting is with the second N in the for loop. When N=0, "print type(tempArr)" returns: type 'numpy.ndarray', which I would expect. When N=1, however, type(tempArr) returns type 'NoneType', which I can't understand. The entire script then fails at line "totalArr += tempArr" because you can't add NoneType to a numpy int array.

What's weird is that no matter how I loop through, the "second" N always fails. For example, "for N in [1, 0, 2]:" returns an array for N=1 and N=2 but returns NoneType when N=0 (even though N=0 didn't fail when going from 0 to Nlayers (aka 0, 1, 2)) and "for N in [2, 1, 0]:" works for N=2 and N=0 but not N=1.

I'm pretty sure there's nothing wrong with the vrt itself, as 'gdalinfo %s' % VRTstack' returns accurate information on all three bands.

Also worth mentioning is that this script worked with one set of images (call it S1), but not another (S2). The only real difference between each set is that S1 has a coarser resolution, and its rasters are 32-bit signed type, while S2's raster type is 8-bit unsigned. I thought that maybe it had something to do with the fact that I was trying to read in unit8 rasters and add it to a uint16, so I tried making totalArray type uint 8 as well, but it still failed. It later occurred to me that this wouldn't change anything, because tempArr is NoneType before it even gets added to totalArray.


When using 5 images instead of just 3, NoneType is always/only returned for the 2nd and 4th N in the loop. That is, when looping through N=0 to N=4, numpy arrays are returned for N=0,2,4, while NoneType is returned for N=1 and N=3. Similarily, when the loop is "for N in [1, 0, 3, 2, 4]," arrays are correctly returned for N=1,3,4 but NoneType is returned for N=0 and N=2. Bizarre.

  • Try adding gdal.UseExceptions() at the top to see if there is more info. I would read all the bands and sum them along the band dimension with one line: vrtDS.ReadAsArray().sum(axis=0) – Mike T Feb 26 '15 at 5:07
  • Thanks Mike. I have gdal.UseExceptions() already and there is no info. The ReadAsArray() doesn't technically fail though, just returns NoneType. I will try your suggestion- I always love a good one-liner to replace a loop, so thank you. – user20408 Feb 26 '15 at 15:31
  • Mike T, I tried your suggestion and am getting "ValueError: total size of new array must be unchanged" at vrtDS.ReadAsArray(). A quick internet search returns nothing useful. Any thoughts? – user20408 Feb 26 '15 at 15:44
  • @MikeT, I don't know if you saw my comment to you since I didn't tag you. If you did see it, sorry for the extra notification. – user20408 Feb 26 '15 at 20:36
  • My guess is there is an issue with the VRT. Do each band overlap exactly and have the same size? If not, you may need to use gdal_merge.py instead. – Mike T Feb 26 '15 at 22:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.