Community

Comparing DEMs with RMSE

1236
1
02-14-2020 09:42 PM
ElsaSartor
by
New Contributor
‎02-14-2020 09:42 PM

Hi everyone! I am attempting to compare 5 DEMs I interpolated with different methods from a point file of elevations with another DEM that will act as the true elevation. I know I need to compute the RMSE to compare each of my interpolated DEMs with the other DEM but I have no idea how to begin going about this in ArcMap. I am required to produce a tool in Modelbuilder to do this comparison as well. I am completely lost on how to start, so any tips would be greatly appreciated!!

I am using ArcMap 10.7.1. I have all my rasters projected into the same spatial reference (NAD 1983 UTM Zone 17) and they are all clipped to the same geographical extent. I have stored the 6 DEMs mentioned above in a file geodatabase.

Thank you in advance for any help that anyone can provide!

0 Kudos
1 Reply
DavidPike
by MVP Frequent Contributor
MVP Frequent Contributor
‎02-15-2020 01:42 AM

Raster calculator the difference between the true DEM and your interpolated one, square it, the run something like zonal stats as table to grab the mean, then take the root.

here's some code. I've not tested it. may have to save the raster objects as interim data. let me know if it doesn't work and i'll amend the code.

import arcpy

#path to your original accurate DEM - make sure path is enclosed by r' '
original_dem = r'c:\\myproject\mygeodatabase\originalDEM'

#change the below paths to your interpolated DEM paths 
# - make sure path is enclosed by r' '
#interpolated DEM 1
interp_dem1 = r'c:\\myproject\mygeodatabase\interpDEM1'
#interpolated DEM 2
interp_dem2 = r'c:\\myproject\mygeodatabase\interpDEM2'
#interpolated DEM 3
interp_dem3 = r'c:\\myproject\mygeodatabase\interpDEM3'
#interpolated DEM 4
interp_dem4 = r'c:\\myproject\mygeodatabase\interpDEM4'
#interpolated DEM 5
interp_dem5 = r'c:\\myproject\mygeodatabase\interpDEM5'

#make a list of the interpolated DEM paths
interp_dem_list = [interp_dem1, interp_dem2, interp_dem3, interp_dem4, interp_dem5]

#turn the original DEM into a raster object,
#iterate through the list, turn paths into raster objects to perform
#raster algebra on them, output raster cell values will be
#(interp_dem - original_dem) ^ 2
original_dem_rasobj = arcpy.Raster(original_dem)
#set a counter to label result with iteration number
iteration_counter = 1
for interp_dem_path in interp_dem_list:
    ras_object = arcpy.Raster(interp_dem_path)
    ras_diff = (ras_object - original_dem_rasobj) ** 2
    #get the mean of this raster
    ras_stats = arcpy.CalculateStatistics_management(ras_diff)
    ras_prop = arcpy.GetRasterProperties_management(ras_stats, "MEAN")
    ras_mean = ras_prop.getOutput(0)
    rmse = float(ras_mean) ** 0.5
    #print out the rmse in the python window
    print("RMSE of Interpolated DEM Raster " + str(iteration_counter) + " = " + str(rmse))
    #add 1 to the iteration counter
    iteration_counter += 1

print("Processing Complete")

‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍
0 Kudos