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()

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

Answer 1

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:

enter image description here

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

enter image description here

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

Answer 2

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)