Getting IndexError while classifying .tif using GDAL python

Question

In the following code (derived from code written by dmh126) I want to classify .tif file based on it's pixel values after converting it into an array. I think there is some mistake in for loop which is giving me an error "IndexError: index 101 is out of bounds for axis 0 with size 101". I am posting my code here.

from osgeo import gdal
import numpy as np
driver = gdal.GetDriverByName( 'GTiff')
file = gdal.Open( 'Gadchiroli.tif')
band = file.GetRasterBand(1)
lista = band.ReadAsArray(0,0,file.RasterXSize,file.RasterYSize).astype(np.int)
print(lista)
print(range(file.RasterXSize))
print(range(file.RasterYSize))

# reclassification
for j in range(file.RasterXSize-1):
    s=lista[j]
    print(s)
    i=0
    #print(j,end='')
    for i in  range(file.RasterYSize-1):
        print(i,j,end='')
        if s[i] < 280:
           print(s[i])
           lista[i,j] = 1
        elif 280 < s[i]< 290:          
            lista[i,j] = 2
        elif 290 < s[i] < 300:
            lista[i,j] = 3
        elif 300 < s[i] < 310:
            lista[i,j] = 4
        elif 310 < s[i] < 320:
          #print(s[i])
            lista[i,j] = 5
        elif 320 < s[i] < 330:
            lista[i,j] = 6
       # print(s[i])
        else:
            lista[i,j] = 7

# create new file
file2 = driver.Create( 'class.tif', file.RasterXSize , file.RasterYSize , 1)
file2.GetRasterBand(1).WriteArray(lista)

# spatial ref system
proj = file.GetProjection()
georef = file.GetGeoTransform()
file2.SetProjection(proj)
file2.SetGeoTransform(georef)
file2.FlushCache()

It doesn't work because you miss datatype parameter in 'Create' method and close file2. Please, see my answer. — xunilk, Jan 06 '18 at 13:18
You have failed to credit the original author of this code as required by the Creative Commons Attribution Share Alike license — user2856, Jan 07 '18 at 01:02

score 1 · Answer 1 · answered Jan 06 '18 at 13:16

It doesn't work because you need to write this line:

file2 = driver.Create( 'class.tif', file.RasterXSize , file.RasterYSize , 1)

as (you miss data type):

file2 = driver.Create( 'pyqgis_data/class.tif', 
                       file.RasterXSize , 
                       file.RasterYSize , 
                       1,
                       gdal.GDT_Int32)

and you also miss close file2.

Following condensed version of your code (without any prints and limits of loops fixed) works:

from osgeo import gdal
import numpy as np
driver = gdal.GetDriverByName( 'GTiff')
file = gdal.Open( 'pyqgis_data/Gadchiroli.tif')
band = file.GetRasterBand(1)
lista = band.ReadAsArray(0,0,file.RasterXSize,file.RasterYSize).astype(np.int)

# reclassification
for j in range(file.RasterXSize):

    s = lista[j]
    i = 0
    for i in  range(file.RasterYSize):
        if s[i] < 280:
           lista[i,j] = 1
        elif 280 < s[i]< 290:          
            lista[i,j] = 2
        elif 290 < s[i] < 300:
            lista[i,j] = 3
        elif 300 < s[i] < 310:
            lista[i,j] = 4
        elif 310 < s[i] < 320:
            lista[i,j] = 5
        elif 320 < s[i] < 330:
            lista[i,j] = 6

        else:
            lista[i,j] = 7

print lista

file2 = driver.Create( 'pyqgis_data/class.tif', 
                       file.RasterXSize , 
                       file.RasterYSize , 
                       1,
                       gdal.GDT_Int32)

file2.GetRasterBand(1).WriteArray(lista)

proj = file.GetProjection()
georef = file.GetGeoTransform()
file2.SetProjection(proj)
file2.SetGeoTransform(georef)
file2.FlushCache()

file2 = None

I tried it out with my own version of Gadchiroli.tif: a random raster of 400x400 with values between 1 and 400. It looks like:

After running code at Python Console, I got this one:

However, you should review classification code because it has serious issues; but it is a subject for another question.

user2856 · Answer 2 · 2018-01-07T05:35:56.847

I would try numpy array operations rather than looping. Much quicker. See below:

from osgeo import gdal
import numpy as np

x=np.random.randint(350, size=(400,500))

def a():
    y = np.zeros(x.shape, dtype=x.dtype) + 7
    y[x < 280]  = 1
    y[(280 < x) & (x < 290)] = 2
    y[(290 < x) & (x < 300)] = 3
    y[(300 < x) & (x < 310)] = 4
    y[(310 < x) & (x < 320)] = 5
    y[(320 < x) & (x < 330)] = 6

def b():
    y = np.zeros(x.shape, dtype=x.dtype) + 7  
    for i in range(x.shape[0]):
        for j in range(x.shape[1]):
            if x[i,j] < 280:
               y[i,j] = 1
            elif 280 < x[i,j]< 290:          
                y[i,j] = 2
            elif 290 < x[i,j] < 300:
                y[i,j] = 3
            elif 300 < x[i,j] < 310:
                y[i,j] = 4
            elif 310 < x[i,j] < 320:
                y[i,j] = 5
            elif 320 < x[i,j] < 330:
                y[i,j] = 6


%timeit a()
    2.82 ms ± 16.9 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

%timeit b()
    159 ms ± 2.63 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

Getting IndexError while classifying .tif using GDAL python

2 Answers2