topic Relatively complex mathematical process on raster grid - beginner in Python Questions

Relatively complex mathematical process on raster grid - beginner

MichaelHammond — Sat, 26 Oct 2013 13:41:23 GMT

Hello, I'm a beginner in ArcGIS / Python and trying to perform a calculation on a raster grid.

I have a raster grid and I want to do some calculations on it using some relatively complex (more so than add, multiply, etc). I've worked out the python code I'd need to do were it not a raster grid.

a = 0.5*(1+math.erf((RASTER)/math.sqrt(2*2*2))) where DATA is the raster I want to put in

Is this something I am better doing with a Python Script, one that imports math? Would anyone know a good point to start please?

Thanks

Re: Relatively complex mathematical process on raster grid - beginner

Luke_Pinner — Sun, 12 Dec 2021 00:40:15 GMT

You can't pass a raster to the math.erf function. You'll get a "TypeError: a float is required" exception.

Here's some python code to run your expression:
Note - I take no responsibility for the correctness of the ERF implementation! I just copied it off teh interwebz and modified it to work with raster objects...

import arcview,arcpy,math
from arcpy.sa import *

def erf(x):
    ''' Error function modified for rasters from the pure python implementation
        posted on StackExchange [1] of the algorithm from the
        Handbook of Mathematical Functions [2]
        [1] http://stackoverflow.com/a/457805/737471
        [2] http://people.math.sfu.ca/~cbm/aands/frameindex.htm
    '''

    # save the sign of x
    #sign = 1 if x >= 0 else -1 # This doesn't work with a Raster object
                                #   ValueError: The truth value of a raster is ambiguous.
                                #   Invalid use of raster with Boolean operator or function.
    sign = (x>=0) + ((x<0)*-1)
    x = Abs(x)

    # constants
    a1 =  0.254829592
    a2 = -0.284496736
    a3 =  1.421413741
    a4 = -1.453152027
    a5 =  1.061405429
    p  =  0.3275911

    # A&S formula 7.1.26
    t = 1.0/(1.0 + p*x)
    y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*Exp(-x*x)
    return sign*y

if __name__=='__main__':
    input_raster='path/to/input/raster'
    output_raster='path/to/output/raster'

    arcpy.CheckOutExtension('spatial')

    inras=Raster(input_raster)
    outras = 0.5*(1+erf((inras)/math.sqrt(2*2*2)))
    outras.save(output_raster)

Re: Relatively complex mathematical process on raster grid - beginner

DouglasSands — Sun, 12 Dec 2021 00:40:18 GMT

I think you can simplify the denominator. The square root of 2 cubed:

math.sqrt(2*2*2)

will always be a constant (2.828...).