[GIS] Raster reclassify using python, gdal and numpy

gdalnumpypython-2.7reclassify

I would like to reclassify a raster file from a raster with 10 classes to a raster with 8 classes using Python, GDAL and/or numpy. The classes are represented as integers. I have tried following the steps from this post Reclassifying rasters using GDAL and Python? , the numpy.equal doc and also gdal_calc doc. However, to no avail.

The raster file to be reclassified has integer values ranging from 0 to 11 and also include values 100 and 255. The following show the reclass (from value : to value):

nodata : 4,
0 : 4,
1 : 1,
2 : 2,
3 : 3,
4 : 3,
5 : 4,
6 : 5,
7 : 5,
8 : 6,
9 : 7,
10 : 8,
100 : nodata,
255 : nodata,

What I have been able to do is select the raster file to be reclassified using tkinter.FileDialog and get the raster info such as geotransform, and pixel size with reclass = gdal.Open(raster, GA_ReadOnly).

How do I go about solving the above.

It might be worth mentioning that the rasters to be reclassified can be fairly large in some cases (500mb to 5gb).

Best Answer

Instead of doing the reclassification as a double for loop described by dmh126, do it using np.where:

# reclassification    
lista[np.where( lista < 200 )] = 1
lista[np.where((200 < lista) & (lista < 400)) ] = 2
lista[np.where((400 < lista) & (lista < 600)) ] = 3
lista[np.where((600 < lista) & (lista < 800)) ] = 4
lista[np.where( lista > 800 )] = 5

On an array of 6163 by 3537 pixels (41.6mb) the classification is done in 1.59 seconds, where it takes 12min 41s using the double for loop. In total just a speedup of 478x.

Bottomline, never use a double for loop using numpy

Related Question