1

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?

Community
  • 1
Ralu Bur
  • 573
  • 5
  • 24
  • 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'. – Mikkel Lydholm Rasmussen Jul 08 '16 at 11:40
  • @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. – Ralu Bur Jul 08 '16 at 12:10

1 Answers1

4

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
Juan
  • 353
  • 2
  • 6