# Estimating slope by fitting a quadratic surface (Evans 1980) in ArcGIS Pro

I am needing to calculate slope using the quadratic surface method of Evans (An integrated system of terrain analysis and slope mapping, 1980), but not sure how to implement it in ArcGIS Pro. A GRASS version of the code is located in this thesis, and the source code is included in GRASS as r.param.scale, but I don't really speak GRASS or C.

In ArcGIS, slope is estimated as explained here. Essentially the equation is the sqrt(rise^2 + run^2). This is different from how Evans (1980) estimated slope.

Evans (1980) estimated slope from a fitted surface, with the formula

z = ax^2 + by^2 + cyx + dx + ey + f, with a,b,c,d,e,f as constants.

slope is sqrt(d^2 + e^2), while curvature is 2a + 2b.

Does anyone have any knowledge as to how to implement this equation in ArcGIS Pro?

I apologize for a lack of a reproducible example, but if someone points me in the in the right direction, I will definitely add a reproducible example.

• I really do not understand your question "how this is implemented in ArcGIS". The linked documentation is fairly clear and the associated "slope" and "curvature" tools implement the method. pro.arcgis.com/en/pro-app/help/data/imagery/slope-function.htm – Jeffrey Evans Mar 12 at 19:43
• I have edited my question to hopefully make it more explicit that ArcGIS actually uses a different estimation of slope than the Evans (1980) paper. The curvature estimate in ArcGIS fits a 9 term polynomial while slope is derived using a different methodology. I would like to implement the methodology of Evans (1980), but just not sure how to go about doing it. Plus, I can't find the paper. – user44796 Mar 12 at 20:30
• Do you know how to use python script? – FelixIP Mar 13 at 8:23
• Yes, that is what I am assuming how I would have to code it up. – user44796 Mar 13 at 13:01
• I agree, and that will likely be my first approach. You are talking about a purely python approach if I am converting to a numpy Array. I think that I should be able to handle it in ArcGIS with raster calculator, just not sure how to get started coding up the underlying algorithm. – user44796 Mar 13 at 14:24

The script below finds 6 coefficients of the model using least squares technique. I tested it using surface computed as

``````Z=0.1*D + 0.05*D*D
``````

where D is euclidean distance. Slope (S) for that surface is known:

``````S=0.1+0.1*D
``````

Real life output example:

``````import arcpy
import numpy as np

##    find 6 parameters of surface equation
def createArrays(validList,zList):
aList=[]
for i in validList:
XX=X[i];YY=Y[i]
aList.append([XX*XX,YY*YY,XX*YY,XX,YY,1])
matrix=[]
yS=[]
for column in range(6):
bList=[0.]*6
igrek=0.
for row in range(len(validList)):
aRow=aList[row]
m=aRow[column]
for n,item in enumerate(aRow):
bList[n]+=item*m
igrek+=zList[row]*m
matrix.append(bList)
yS.append(igrek)
matrix=np.array(matrix);yS=np.array(yS)
matrix=matrix.T
coeffs=np.linalg.solve(matrix,yS)
return coeffs

##    cell coordinates relative to one in question
X=[-1,0,1,-1,0,1,-1,0,1]
Y=[1,1,1,0,0,0,-1,-1,-1]
##    get dem raster and compute dimensions
cSize=float(zRaster.meanCellHeight)
zArray = arcpy.RasterToNumPyArray(zRaster,"","","",-999)
nRows,nCols=zArray.shape
cellsTotal=nCols*nRows
blankArray=arcpy.RasterToNumPyArray(zRaster,"","","",-999)
arcpy.SetProgressor("step", "", 0, cellsTotal,1)
##    iterate through Z raster
for nRow in range(nRows):
for nCol in range(nCols):
arcpy.SetProgressorPosition()
##    get cell neighbours with data and fill lists for least squares
n2pick=-1
vList=[]
zS=[]
for iY in range(-1,2):
for iX in range(-1,2):
n2pick+=1
nR=nRow+iY
if nR not in range(nRows):continue
nC=nCol+iX
if nC not in range(nCols):continue
theZ=zArray[nR,nC]
if theZ==-999:continue
vList.append(n2pick)
zS.append(theZ)
if len(vList)<7:
blankArray[nRow,nCol]=None
else:
A,B,C,D,E,F=createArrays(vList,zS)
blankArray[nRow,nCol]=math.hypot(D,E)/cSize
##  convert results to grid
d=arcpy.Describe(zRaster)
origin=d.extent.lowerLeft
myRaster = arcpy.NumPyArrayToRaster(blankArray,origin,cSize,cSize)
outRaster = arcpy.GetParameterAsText(0)
myRaster.save(outRaster)
``````
• That worked perfectly. Just to be clear, the output is rise/run, right? BTW, I compared the results to r.param.scale in GRASS and it seems like your implementation is providing more sensitive results. r.param.scale was classifying areas I know are flat as having quite a bit of slope. – user44796 Mar 22 at 14:32
• Yes it does compute slope for that type of surface approximation. – FelixIP Mar 22 at 19:49