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.

[DarrenWiens2]

Have you considered flipping to Numpy, reclassifying, and then back again?

RasterToNumPyArray—ArcPy Functions | ArcGIS for Desktop

NumPyArrayToRaster—ArcPy Functions | ArcGIS for Desktop

edit: I assume you've considered and passed on using the Reclassify geoprocessing tool.

[DarrenWiens2]

Here is approximately how you can do this in numpy (need to do some other tedious stuff to position final raster correctly):

>>> import numpy   
... reclass_map = [[10,999],[20,1000]] # reclass list: [[old1,new1],[old2,new2],...]
... arr = arcpy.RasterToNumPyArray('idw.tif',nodata_to_value=0) # raster to  numpy
... for i in reclass_map: # loop through class list
...    temp = numpy.less(arr,i[0]) # if less than value
...    numpy.putmask(arr,temp,i[1]) # set to new value
... new_ras = arcpy.NumPyArrayToRaster(arr) # flip back to raster

[DanPatterson_Retired]

It seems the original post got moved or something... and my mention deleted, so here it is again

Conceptually and overtly just convert to array first, then zz back to raster

import numpy as np
a = np.arange(60).reshape((6,10))
zz = np.zeros_like(a)
bins = [0,5,10,15,20,25,30,60,100]
newclass = [1, 2, 3, 4, 5, 6, 7, 8]
noob = zip(bins[:-1],bins[1:],newclass)
qs =[]
for rc in noob:
    q1 = (a >= rc[0])
    q2 = (a < rc[1])
    z = np.where(q1 & q2, rc[2],0)
    zz = zz + z
print("old\n{}\nnew\n{}".format(a,zz))


old
[[ 0  1  2  3  4  5  6  7  8  9]
[10 11 12 13 14 15 16 17 18 19]
[20 21 22 23 24 25 26 27 28 29]
[30 31 32 33 34 35 36 37 38 39]
[40 41 42 43 44 45 46 47 48 49]
[50 51 52 53 54 55 56 57 58 59]]
new
[[1 1 1 1 1 2 2 2 2 2]
[3 3 3 3 3 4 4 4 4 4]
[5 5 5 5 5 6 6 6 6 6]
[7 7 7 7 7 7 7 7 7 7]
[7 7 7 7 7 7 7 7 7 7]
[7 7 7 7 7 7 7 7 7 7]]
old
[[ 0  1  2  3  4  5  6  7  8  9]
[10 11 12 13 14 15 16 17 18 19]
[20 21 22 23 24 25 26 27 28 29]
[30 31 32 33 34 35 36 37 38 39]
[40 41 42 43 44 45 46 47 48 49]
[50 51 52 53 54 55 56 57 58 59]]
new
[[1 1 1 1 1 2 2 2 2 2]
[3 3 3 3 3 4 4 4 4 4]
[5 5 5 5 5 6 6 6 6 6]
[7 7 7 7 7 7 7 7 7 7]
[7 7 7 7 7 7 7 7 7 7]
[7 7 7 7 7 7 7 7 7 7]]

[DarrenWiens2]

Your mention was just me (cheekily) alerting you that there was an opportunity to use numpy, but after I provided one such solution, it didn't seem necessary.

[DanPatterson_Retired]

naw... it just said "mention removed" on the mail link I had to post it twice, since you were probably editing yours when I posted and My original didn't keep

[PeterWilson]

Hi Dan

The following looks to be what I was looking for. I've been swamped the last few days. Will come back to the following over the weekend and post my code using numpy.

[PeterWilson]

Hi Dan

Your sample code worked perfectly, with some minor changes.

'''
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:\Python\Testing\Numpy_Reclassify\SAA_Raster\tbr")
output_folder = r"E:\Python\Testing\Numpy_Reclassify\SAA_Reclassified"


# 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, 5, 10, 15, 20, 25, 30, 60]
    new_bins = [5, 10, 15, 20, 25, 30, 60]
    new_classes = zip(bins[:-1], bins[1:], new_bins)
    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)
    Int(reclass_raster).save(os.path.join(output_folder, "rec"))


reclass_raster(input_raster, output_folder)


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

Regards

[DanPatterson_Retired]

looks good Peter

[SteveLynch]

Peter Wilson​ why did you not consider any of Spatial Analyst's reclass tools? They now support muti-core processing

-Steve