1

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.

Example of classifying quantile breaks in ArcMap

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)
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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.