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

Question

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?

Answer 1

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