Reclassifying raster using numpy

Question

I have more than 300 NDVI, Soil Salinity and NDMI Calculated Rasters I want to reclassify into categories based on percentages of values like this function:

def classify_salinity(ras):
 max_val = np.nanmax(ras)
 min_val = np.nanmin(ras)
 #scaling between 0-100%
 scaled = ((ras - min_val) / (max_val - min_val)) * 100
 """
 reclassifying according to critera
 Criteria for Salinity:
 {>80% == 1, 50-80% == 2, 20-50% == 3, <20% == 4}
 """
 # Float to integer and copy array
 arr = scaled
 arr[(arr < 20)] = 4
 arr[(arr >= 80)] = 1
 arr[(arr >= 50) & (arr < 80)] = 2
 arr[(arr >= 20) & (arr < 50)] = 3
 # arr[(arr < 0)] = 0
 return scaled

The issue is that this function executes statements line by line on the array which can introduce complications.

Is there a way to do this all in one statement using numpy?

• Kadir Şahbaz: What do you mean by "one statement"? One line of code?
• Hornbydd: Some ideas here.

Answer 1

You can iterate all over array, if that helps?

def reclassify(arr, rows, cols):
 for r in range(rows):
 for c in range(cols):
 temp_value = arr[r][c]
 if temp_value < 20:
 temp_value = 4
 elif temp_value >=80:
 temp_value = 1
 
 elif (temp_value >=50) and (temp_value <80):
 temp_value = 2
 
 elif (temp_value >=20) and (temp_value <50):
 temp_value = 3
 
 else:
 pass
 
 arr[r][c] = temp_value
 
 return arr

Answer 2

It works in my test

 import numpy as np
 from osgeo import gdal
 
 driver = gdal.GetDriverByName('GTiff')
 file = gdal.Open('example.tif')
 band = file.GetRasterBand(1)
 ras = band.ReadAsArray()
 max_val = np.nanmax(ras)
 min_val = np.nanmin(ras)
 
 #scaling between 0-100%
 lista = ((ras - min_val) / (max_val - min_val)) * 100
 
 # reclassification
 list_dest = lista.copy()
 
 list_dest[np.where(lista >= 80)] = 1
 list_dest[np.where((lista >= 50) & (lista < 80)) ] = 2
 list_dest[np.where((lista >= 20) & (lista < 50)) ] = 3
 list_dest[np.where(lista < 20)] = 4
 # create new file
 file2 = driver.Create( 'example_reclassed.tif', file.RasterXSize , file.RasterYSize , 1)
 file2.GetRasterBand(1).WriteArray(list_dest)
 
 # spatial ref system
 proj = file.GetProjection()
 georef = file.GetGeoTransform()
 file2.SetProjection(proj)
 file2.SetGeoTransform(georef)
 file2.FlushCache()