demo_reclassify_cases.py

import pygeoprocessing
import numpy
from osgeo import gdal, osr

srs = osr.SpatialReference()
srs.ImportFromEPSG(32731)  # WGS84/UTM zone 31s
projection_wkt = srs.ExportToWkt()
arr = numpy.array([[0, 1, -1]], dtype=numpy.int16)
base_nodata = 0
target_datatype = gdal.GDT_Int16
in_path = '/Users/emily/Documents/reclassify_raster_input.tif'
out_path = '/Users/emily/Documents/reclassify_raster_output.tif'
pygeoprocessing.numpy_array_to_raster(
    arr, base_nodata, (1, 1), (0, 0), projection_wkt, in_path)


######################################################################
# Case 1: nodata not in value_map, target_nodata same as input nodata
value_map = {
    1: 5,
    -1: 6
}
target_nodata = base_nodata

pygeoprocessing.reclassify_raster(
    (in_path, 1),
    value_map,
    out_path,
    target_datatype,
    target_nodata)

# nodata is preserved
expected_out = numpy.array([[0, 5, 6]])
out_nodata = pygeoprocessing.get_raster_info(out_path)['nodata'][0]
out = pygeoprocessing.raster_to_numpy_array(out_path)
print('\n\nCase 1')
print('Expected values:', expected_out)
print('Actual values:', out)
print('Expected nodata value:', target_nodata)
print('Actual nodata value:', out_nodata)


############################################################################
# Case 2: nodata not in value_map, target_nodata different from input nodata
value_map = {
    1: 5,
    -1: 6
}
target_nodata = 10

pygeoprocessing.reclassify_raster(
    (in_path, 1),
    value_map,
    out_path,
    target_datatype,
    target_nodata)

# nodata pixels are not reclassified
# the output raster's nodata value is changed
expected_out = numpy.array([[0, 5, 6]])
out_nodata = pygeoprocessing.get_raster_info(out_path)['nodata'][0]
out = pygeoprocessing.raster_to_numpy_array(out_path)
print('\n\nCase 2')
print('Expected values:', expected_out)
print('Actual values:', out)
print('Expected nodata value:', target_nodata)
print('Actual nodata value:', out_nodata)

###################################################################
# Case 3: nodata is in value_map, target_nodata is same
value_map = {
    0: 10,
    1: 5,
    -1: 6
}
target_nodata = 10

pygeoprocessing.reclassify_raster(
    (in_path, 1),
    value_map,
    out_path,
    target_datatype,
    target_nodata)

# nodata is reclassified to `target_nodata`
# the output raster's nodata value is `target_nodata`
expected_out = numpy.array([[10, 5, 6]])
out_nodata = pygeoprocessing.get_raster_info(out_path)['nodata'][0]
out = pygeoprocessing.raster_to_numpy_array(out_path)
print('\n\nCase 3')
print('Expected values:', expected_out)
print('Actual values:', out)
print('Expected nodata value:', target_nodata)
print('Actual nodata value:', out_nodata)

##################################################################
# Case 4: nodata is in value_map, target_nodata is different
value_map = {
    0: 10,
    1: 5,
    -1: 6
}
target_nodata = 0

pygeoprocessing.reclassify_raster(
    (in_path, 1),
    value_map,
    out_path,
    target_datatype,
    target_nodata)

# nodata pixels are reclassified to a valid value
expected_out = numpy.array([[10, 5, 6]])
out_nodata = pygeoprocessing.get_raster_info(out_path)['nodata'][0]
out = pygeoprocessing.raster_to_numpy_array(out_path)
print('\n\nCase 4')
print('Expected values:', expected_out)
print('Actual values:', out)
print('Expected nodata value:', target_nodata)
print('Actual nodata value:', out_nodata)