topic Re: Generate a Function for Map Algebra Expression in Python Questions

Generate a Function for Map Algebra Expression

PeterWilson — Fri, 11 Mar 2016 18:42:19 GMT

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.

Re: Generate a Function for Map Algebra Expression

DarrenWiens2 — Fri, 11 Mar 2016 18:47:36 GMT

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.

Re: Generate a Function for Map Algebra Expression

DarrenWiens2 — Sat, 11 Dec 2021 21:49:44 GMT

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

Re: Generate a Function for Map Algebra Expression

DanPatterson_Retired — Sat, 11 Dec 2021 21:49:46 GMT

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

Re: Generate a Function for Map Algebra Expression

DarrenWiens2 — Fri, 11 Mar 2016 21:33:19 GMT

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.

Re: Generate a Function for Map Algebra Expression

DanPatterson_Retired — Fri, 11 Mar 2016 21:38:38 GMT

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

Re: Generate a Function for Map Algebra Expression

PeterWilson — Thu, 17 Mar 2016 08:11:38 GMT

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.

Re: Generate a Function for Map Algebra Expression

PeterWilson — Sat, 11 Dec 2021 21:49:49 GMT

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

Re: Generate a Function for Map Algebra Expression

DanPatterson_Retired — Mon, 28 Mar 2016 21:24:39 GMT

looks good Peter

Re: Generate a Function for Map Algebra Expression

SteveLynch — Mon, 28 Mar 2016 22:20:49 GMT

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

-Steve

Re: Generate a Function for Map Algebra Expression

DanPatterson_Retired — Mon, 28 Mar 2016 22:26:34 GMT

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

Re: Generate a Function for Map Algebra Expression

PeterWilson — Tue, 29 Mar 2016 08:24:31 GMT

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.

Re: Generate a Function for Map Algebra Expression

PeterWilson — Sat, 11 Dec 2021 21:49:52 GMT

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

Re: Generate a Function for Map Algebra Expression

DanPatterson_Retired — Wed, 30 Mar 2016 11:26:49 GMT

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

Re: Generate a Function for Map Algebra Expression

curtvprice — Sat, 11 Dec 2021 21:49:55 GMT

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

Re: Generate a Function for Map Algebra Expression

DanPatterson_Retired — Sat, 09 Apr 2016 05:33:15 GMT

missing the import lines of course

rows and columns of raster?

assuming it is float of course...

just curious

Re: Generate a Function for Map Algebra Expression

curtvprice — Sat, 09 Apr 2016 05:39:57 GMT

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

Re: Generate a Function for Map Algebra Expression

DanPatterson_Retired — Sat, 09 Apr 2016 05:47:29 GMT

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

Re: Generate a Function for Map Algebra Expression

curtvprice — Sat, 09 Apr 2016 05:55:05 GMT

It's not a view but a "temporary raster" in arcpy map algebra... out_raster will be written to a real raster as soon as it's needed in a later map algebra expression. Including if you save it with

out_raster.save("outraster.tif")

In ArcMap at the python prompt, however, it gets calculated and written right away to a temp grid in the scratch folder so it can be added to the map (if addOutputsToMap is set).

The dimensions of the output are determined by the currently active environment settings.

Re: Generate a Function for Map Algebra Expression

DanPatterson_Retired — Sat, 09 Apr 2016 06:44:43 GMT

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

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

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