Classifying quantiles in raster dataset with ArcPy?

Question

I am trying to write a Python script which uses the reclassify function in ArcPy to create 9 classes in a raster file. I want to use the quantile method to create the breaks, but so far, I can't seem to find a function that can create quantile breaks in the raster data.

Essentially, I want to perform the above tool with ArcPy. Currently, I have the below in my code, but I do not know how to calculate the remap variable. It cannot be hardcoded as this function should work for all raster inputs. Also, the only modules I have installed are ArcPy and NumPy. I am not able to install any other modules because this is a school project.

tool_input = in_file
reclass_field = "Value"
remap = ???
missing_values = "NODATA"
tool_output = out_files

Reclassify(tool_input, reclass_field, remap, missing_values).save(tool_output)

Answer 1

Here's a function that will reclassify a raster by quantiles.

I included thorough comments, but some other things to note:

1. numpy.quantile is relatively new, so most ArcGIS releases won't have that function in their numpy versions. Instead, I used numpy.percentile--well, actually, I had to use numpy.nanpercentile to ignore null values.
2. This function will work with integer and float rasters, as well as with on-disk rasters and in-memory arcpy.Raster objects.
3. The "quantiles" argument to the function should simply be an integer indicating the number of output classes you want.

Here are some helpful links for reference:

numpy.nanpercentile

raster to numpy array

remap range

dealing with integer array

import numpy
import arcpy
arcpy.CheckOutExtension('SPATIAL')

def reclassify_by_quantiles(raster, quantiles):
    desc = arcpy.Describe(raster)

    # Get the quantile break points
    percentiles = list()
    for i in range(1, quantiles):
        percentiles.append(i * (100.0/quantiles))

    # Ensure that the raster is a raster object, in order to identify the minimum and
    #   maximum values
    if type(raster) != arcpy.Raster:
        raster = arcpy.Raster(raster)

    value_minimum = raster.minimum
    value_maximum = raster.maximum

    # Identify a value that does not occur in the raster
    null_value = value_minimum - 1

    # If the raster is not an integer type, then there is no need to futz with NaNs
    if desc.pixelType.startswith('F'):
        arr = arcpy.RasterToNumPyArray(raster, nodata_to_value=numpy.NaN)
    # Since integer arrays can't contain NaNs, you must do wackiness
    else:
        # Convert to an array, setting NoData cells to the unique value
        arr = arcpy.RasterToNumPyArray(raster, nodata_to_value=null_value)
        # Convert the array of integers to an array of floats
        arr = arr.astype('float')
        # Replace the placeholder null value with NaNs
        arr[arr==null_value] = numpy.NaN

    # Compile the quantile breaks
    breakpoints = list(numpy.nanpercentile(arr, percentiles))
    breakpoints.insert(0, value_minimum)
    breakpoints.append(value_maximum)

    # You no longer need the array...though you could do the remaining calculations
    #   on the array instead of one the raster
    del(arr)

    # Map the ranges to class numbers
    remap_table = list()
    for index, breakpoint in enumerate(breakpoints[:-1]):
        remap_table.append([breakpoint, breakpoints[index+1], index+1])
    remap = arcpy.sa.RemapRange(remap_table)

    result = arcpy.sa.Reclassify(raster, 'Value', remap)
    return result

input_raster = r'C:\Path\To\My\Input.gdb\Raster'
# Call the function
result = reclassify_by_quantiles(raster, 9)
# Save the result
result.save(r'C:\Path\To\Save\My\Output.gdb\Raster)