Reclassify raster by defined intervals

import rasterio
import numpy as np
import os

print('All libraries successfully imported!')

All libraries successfully imported!

Set directory

computer_path = '/export/miro/ndeffense/LBRAT2104/'
grp_letter    = 'X'

# Directory for all work files
work_path = f'{computer_path}GROUP_{grp_letter}/WORK/'

input_file = f'{work_path}NDVI/T31UFS_20200417T104021_NDVI.tif'

output_file = f'{input_file[:-4]}_reclassified.tif'

Set parameters

nodata_val = -10000

# User must defined intervals
interval = [-1,-0.5,0,0.5,1]

dtype_out = 'int16'

print(interval)

[-1, -0.5, 0, 0.5, 1]

Reclassify raster by defined intervals

Return the indices of the bins to which each value in input array belongs.

right	order of bins	returned index ``i`` satisfies
False	increasing	bins[i-1] <= x < bins[i]
True	increasing	bins[i-1] < x <= bins[i]
False	decreasing	bins[i-1] > x >= bins[i]
True	decreasing	bins[i-1] >= x > bins[i]

By default, right = False

If values in x are beyond the bounds of bins, 0 or len(bins) is returned as appropriate.

https://numpy.org/doc/stable/reference/generated/numpy.digitize.html

Example

x = np.array([1.2, 10.0, 12.4, 15.5, 20.])
bins = np.array([0, 5, 10, 15, 20])

inds_right_true = np.digitize(x,bins,right=True)
inds_right_false = np.digitize(x,bins,right=False)

print(f'Original continuous values : {x}')
print(f'Interval : {bins}')
print('Reclassified discrete values : ')
print(inds_right_true)
print(inds_right_false)

Original continuous values : [ 1.2 10.  12.4 15.5 20. ]
Interval : [ 0  5 10 15 20]
Reclassified discrete values :
[1 2 3 4 4]
[1 3 3 4 5]

src = rasterio.open(input_file, 'r')
im_arr = src.read()

# Update the dtype in raster metadata (= profile)
profile = src.profile
profile.update(dtype = dtype_out)


# Replace -10000 by np.nan
im_arr[im_arr==nodata_val] = np.nan

# Create a mask with all no data value
mask = np.isnan(im_arr)

# Convert interval into array
bins = np.array(interval)

# Return the indices of the bins to which each value in input array belongs
im_arr_reclass = np.digitize(im_arr, bins, right=False)

# Apply mask on reclassified raster
im_arr_reclass = np.where(mask, nodata_val, im_arr_reclass)

# Change dtype of raster to match the dtype of profile
im_arr_reclass = im_arr_reclass.astype(dtype_out)

print(f'data type = {im_arr_reclass.dtype}')

print(im_arr_reclass)

# Write output file
dst = rasterio.open(output_file, "w", **profile)
dst.write(im_arr_reclass)

src.close()
dst.close()

data type = int16
[[[4 4 4 ... 4 4 4]
  [4 4 4 ... 4 4 3]
  [4 4 4 ... 4 3 3]
  ...
  [4 4 4 ... 3 3 3]
  [4 4 4 ... 3 3 3]
  [4 4 4 ... 3 3 3]]]