0

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

0

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)

Your Answer

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

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