Numpy raster reclassify outputs only 2 out of 6 classes

Question

I have a raster (.tif) which I want to classify (6 classes) using gdal, python and numpy. The cell values range from: MIN = -30.9847 to MAX = 7.3505.

From the first two answers in this question, I've come to the following code:

from osgeo import gdal
import numpy as np
driver = gdal.GetDriverByName('GTiff')
input = gdal.Open('D:/test.tif')
band = input.GetRasterBand(1)
matrix = band.ReadAsArray()
# reclassification
matrix[np.where(matrix < 0)] = 1
matrix[np.where((0 < matrix) & (matrix < 0.3)) ] = 2
matrix[np.where((0.3 < matrix) & (matrix < 0.6)) ] = 3
matrix[np.where((0.6 < matrix) & (matrix < 1 )) ] = 4
matrix[np.where((1 < matrix) & (matrix < 1.5)) ] = 5
matrix[np.where((matrix > 1.5)) ] = 6
# create new file
output = driver.Create('D:/test_class.tif', input.RasterXSize, input.RasterYSize, 1)
output.GetRasterBand(1).WriteArray(matrix)
# spatial ref system
proj = input.GetProjection()
georef = input.GetGeoTransform()
output.SetProjection(proj)
output.SetGeoTransform(georef)
output.FlushCache()

The output raster is created, but it only has 2 classes, value 1 and value 6. Why could that happen?

• Mikkel Lydholm Rasmussen: Are you sure that you actually want the values to be float32? It seems like a horrible waste of space to save what is actually a Byte-format image as a float32. If you really want to do float32, add 'gdal.GDT_Float32' to the end of your 'driver.Create'.
• Ralu Bur: @Mikkel, I thought that discarding the float32 values was the problem, since I couldn't correctly generate the classification. I did what you said and you're right, the values now are float32, but I'm still left with only 2 of 6 classes in the output file. I will edit my question accordingly.

Answer 1

You are overwriting your data 'matrix'. First you reclassify the negative to 1, then you do other steps... but at the end yo are doing :

matrix[np.where((matrix > 1.5)) ] = 6

Which is reclassifying your already changed values to 6.

That is why you have only 1 (negative) and 6 the other values (2,3,4,5) -> 6

I would suggest to copy your matrix like this:

import copy
new_matrix = copy.copy(matrix)
# reclassification
new_matrix[np.where(matrix < 0)] = 1
new_matrix[np.where((0 < matrix) & (matrix < 0.3)) ] = 2
new_matrix[np.where((0.3 < matrix) & (matrix < 0.6)) ] = 3
new_matrix[np.where((0.6 < matrix) & (matrix < 1 )) ] = 4
new_matrix[np.where((1 < matrix) & (matrix < 1.5)) ] = 5
new_matrix[np.where((matrix > 1.5)) ] = 6

• Ralu Bur: I didn't realize that. You are right, now it's working like a charm.
• Juan: Happy to help! cheers