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.– SethinacanMar 14, 2016 at 8:00
-
Of what do you want to compute the RMSE? (Elevations? Locations? 3D coordinates?) What exactly does the shapefile represent?– whuberMar 14, 2016 at 16:39
1 Answer
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:
print "Field already exists"
#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 "Add field."
GBGAddField(TheFeatureClassOfPoints, RMSEField, "DOUBLE")
print "Calculating field."
arcpy.CalculateField_management(TheFeatureClassOfPoints, RMSEField, rmse, "VB")
print "Finished"
except arcpy.ExecuteError:
# Get the tool error messages
msgs = arcpy.GetMessages(2)
# Return tool error messages for use with a script tool
arcpy.AddError(msgs)
# Print tool error messages for use in Python/PythonWin
print msgs
except:
# Get the traceback object
tb = sys.exc_info()[2]
tbinfo = traceback.format_tb(tb)[0]
# Concatenate information together concerning the error into a message string
pymsg = "PYTHON ERRORS:\nTraceback info:\n" + tbinfo + "\nError Info:\n" + str(sys.exc_info()[1])
msgs = "ArcPy ERRORS:\n" + arcpy.GetMessages(2) + "\n"
# Return python error messages for use in script tool or Python Window
arcpy.AddError(pymsg)
arcpy.AddError(msgs)
# Print Python error messages for use in Python / Python Window
print pymsg + "\n"
print msgs