# How to calculate and compare RMSE between DEM using different interpolation methods?

I have data from aerial photo, it was built from stereo photo so I have 3D information that I am using for the DEM. I extracted it into shapefile. I want to get RMSE (Root Mean Square Error) from different interpolation data using ArcGIS. How would I go about solving this problem, for interpolation I plan to use TIN, Kriging and IDW, are these the correct tools? Also do I need more data?

• TIN should have 0 error at the points used to create it. – Sethinacan Mar 14 '16 at 8:00
• Of what do you want to compute the RMSE? (Elevations? Locations? 3D coordinates?) What exactly does the shapefile represent? – whuber Mar 14 '16 at 16:39

Yes, you do need more data, you need about 30 point locations and those point locations need to have known (surveyed) elevation values. If you have those points you could add the elevation values from your interpolated surfaces to the point file as a new attribute using Get Surface Information.

I have included some Python code to calculate the RMSE between the two columns of data, a new column gets added that will store the RMSE between the known elevation values and the interpolated values.

``````  try:
import arcpy, exceptions, math, traceback, sys

def GBGAddField(in_table, field_name, field_type, field_precision = None, field_scale = None, field_length = None, field_alias = None, field_is_nullable = None, field_domain = None):
#Adds a new field to a feature class or table.  If the table already exists, this function will end
#gracefully and not crash your script.  Parameters 1,2 and 3 are required, all else are optional.
fieldList = arcpy.ListFields(in_table)
fieldList2 = []
for field in fieldList:
fieldList2.append(field.name)
if field_name not in fieldList2:
arcpy.AddField_management(in_table, field_name, field_type, field_precision, field_scale, field_length, field_alias, field_is_nullable, field_domain)
else:

#TheFeatureClassOfPoints = arcpy.GetParameterAsText(0)
#InterpolatedValueField = gp.GetParameterAsText(1)
#SurveyedValue = gp.GetParameterAsText(2)
#RMSEField = gp.GetParameterAsText(3)

TheFeatureClassOfPoints = r"SamplePoints"
InterpolatedValueField = r"trend_xp"
SurveyedValue = r"WSE"
RMSEField = r"RMSE_dxp"

SquaredValues = []

print "Start cursor"
rows = arcpy.SearchCursor(TheFeatureClassOfPoints)
for row in rows:
#while row:
i = row.getValue(InterpolatedValueField)
s = row.getValue(SurveyedValue)
SquaredValues.append(pow((i - s), 2))

del rows

n = len(SquaredValues)

SumOfSquares = sum(SquaredValues)

x = SumOfSquares/n

rmse = str(pow(x, .5))
print "rmse =", rmse

print "Calculating field."
arcpy.CalculateField_management(TheFeatureClassOfPoints, RMSEField, rmse, "VB")

print "Finished"

except arcpy.ExecuteError:
# Get the tool error messages
msgs = arcpy.GetMessages(2)
# Print tool error messages for use in Python/PythonWin
print msgs
except:
# Get the traceback object
tb = sys.exc_info()
tbinfo = traceback.format_tb(tb)
# Concatenate information together concerning the error into a message string
pymsg = "PYTHON ERRORS:\nTraceback info:\n" + tbinfo + "\nError Info:\n" + str(sys.exc_info())
msgs = "ArcPy ERRORS:\n" + arcpy.GetMessages(2) + "\n"
# Return python error messages for use in script tool or Python Window