Generate a Function for Map Algebra Expression

3536
20
03-11-2016 10:42 AM
PeterWilson
Regular Contributor

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.

0 Kudos
20 Replies
DarrenWiens2
MVP Honored Contributor

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.

0 Kudos
DarrenWiens2
MVP Honored Contributor

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
0 Kudos
DanPatterson_Retired
MVP Esteemed Contributor

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]]
0 Kudos
DarrenWiens2
MVP Honored Contributor

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.

0 Kudos
DanPatterson_Retired
MVP Esteemed Contributor

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

0 Kudos
PeterWilson
Regular Contributor

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
Regular Contributor

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

0 Kudos
DanPatterson_Retired
MVP Esteemed Contributor

looks good Peter

0 Kudos
SteveLynch
Esri Regular Contributor

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

-Steve

0 Kudos