Generate a Function for Map Algebra Expression

[PeterWilson]

I've written a Python script that reclassifies a raster using a set of conditional statements into 4 classes:

• 0 - 12
• 12 - 24
• 24 - 30
• 30 - 60

rec_ras = Con(costmin <= 12, 12, Con((costmin > 12) & (costmin <= 24), 24, Con((costmin > 24) & (costmin <= 30), 30, Con((costmin > 30) & (costmin <= 60), 60))))

I then needed to change the classes into 7 classes:

• 0 - 5
• 5 - 10
• 10 - 15
• 15 - 20
• 20 - 25
• 25 - 30
• 30 - 60

rec_ras = Con(costmin <= 5, 5, Con((costmin > 5) & (costmin <= 10), 10, Con((costmin > 10) & (costmin <= 15), 15, Con((costmin > 15) & (costmin <= 20), 20, Con((costmin > 20) & (costmin <= 25), 25, Con((costmin > 25) & (costmin <= 30), 30, Con((costmin > 30) & (costmin <= 60), 60)))))))

Is ther a better way of creating a Python function that can create my conditional statement as and if when I need to change the classes before running it. The following is part of  a larger Python Class.

[DanPatterson_Retired]

Steve... I venture is was because of the speed and pure beauty of numpy  

[PeterWilson]

Hi Dan

I found that the new script that I wrote using Numpy has a further advantage in that it supports reclassifying an integer raster to float raster.

'''
Created on Mar 19, 2016

@author: PeterW
'''
# import site-packages and modules
import os
import numpy as np
import arcpy
from arcpy.sa import *

# set arguments
input_raster = arcpy.Raster(r"E:\Projects\2016\G112030\SAA\Rasters\euc_road_clip")
output_folder = r"E:\Projects\2016\G112030\SAA\Rasters"

# set environment settings
arcpy.env.overwriteOutput = True
arcpy.env.outputCoordinateSystem = input_raster

# check out extensions
arcpy.CheckOutExtension("Spatial")


def reclass_raster(input_raster, output_folder):
    lower_left = arcpy.Point(input_raster.extent.XMin,
                            input_raster.extent.YMin)
    cell_size = input_raster.meanCellWidth
    saa_array = arcpy.RasterToNumPyArray(input_raster)
    reclass_array = np.zeros_like(saa_array)
    bins = [0, 10, 25, 50, 100, 13069.5234375]
    new_bins = [1, 1.25, 1.5, 2, 5]
    new_classes = zip(bins[:-1], bins[1:], new_bins)
    print(new_classes)
    for rc in new_classes:
        q1 = (saa_array >= rc[0])
        q2 = (saa_array < rc[1])
        z = np.where(q1 & q2, rc[2], 0)
        reclass_array = reclass_array + z
    reclass_raster = arcpy.NumPyArrayToRaster(reclass_array,
                                              lower_left,
                                              x_cell_size=cell_size,
                                              value_to_nodata=0)
    reclass_raster.save(os.path.join(output_folder, "cost_road"))

reclass_raster(input_raster, output_folder)

# check in extensions
arcpy.CheckInExtension("Spatial")

[DanPatterson_Retired]

Nice...  further reducing steps   and showing how arcpy and numpy get along well.

[PeterWilson]

Hi Steve

I was originally using map algebra Conditional Statement to reclassify my output raster and was looking for a way to build a function for the new classes that I wanted the raster to be reclassified to. Darren and Dan suggested that Numpy would allow me to generate a function to define the classes more cleanly. I must admit I had not considered the reclassify tool, RemapRange would have worked just just as well.

[curtvprice]

You numpy fanboys. I mean really. 

Another thing about spatial analyst tools is they are somewhat more likely to work when you throw a 5GB raster at them. As I have been discovering recently processing big chunks of 10m elevation data!

RemapRange—Spatial Analyst module | ArcGIS for Desktop

Four lines, my friends! (OK, Dan, six lines with the imports)

>>> import arcpy
>>> from arcpy.sa import *
>>> bins = [0, 10, 25, 50, 100, 13069.5234375] 
>>> new_bins = [1, 1.25, 1.5, 2, 5] 
>>> remap = [[bins, bins[k+1], int(new_bins * 100)] for k in range(len(bins) - 1)]
>>> for kk in remap: print(kk)
...
[0, 10, 100]
[10, 25, 125]
[25, 50, 150]
[50, 100, 200]
[100, 13069.5234375, 500]
>>> out_raster =  Reclassify("in_raster", "VALUE", RemapRange(remap)) * 0.01

[DanPatterson_Retired]

missing the import lines of course

rows and columns of raster?

assuming it is float of course...

just curious

[curtvprice]

Okay, I added the import lines for you. This is Python map algebra, I don't need to worry about the raster dimensions!

I really liked Peter's nifty construction of the remap table below, but I needed to scale the data for the float output he was looking for, so I had to go with my somewhat clunkier (but still so pythonesque) list comprehension.

new_classes = zip(bins[:-1], bins[1:], new_bins)  # cool

[DanPatterson_Retired]

Oh... it just changes the 'view' of the raster, it doesn't actually produce a new reclassed one (if that is correct)

[DanPatterson_Retired]

yes, one of several ways to generate it using LCs, but be cautious of zip

When people start using python 3.x + make sure to use  list(zip(....))  or you will get a zip object... which is useful, or not depending upon what you need to do.

bins  = [0, 10, 25, 50, 100, 13079]

new_bins = [1, 1.25, 1.5, 2, 5]

nc = [ (bins[i-1], bins, new_bins[i-1]) for i in range(1,len(bins)) ]  # safe in python 2.x and 3.x

or

zip(bins[:-1],bins[1:],new_bins)  # python 2.x only

or

list(zip(blah)) # in python 3.x

yields

[(0, 10, 1), (10, 25, 1.25), (25, 50, 1.5), (50, 100, 2), (100, 13079, 5)]

[curtvprice]

Sounds like yet another thing to do slightly differently to make your code work in both.

list(zip(blah)) it is.