Comparing DEMs with RMSE

ElsaSartor

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!

DavidPike

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.

<SPAN class="keyword token">import</SPAN> arcpy

<SPAN class="comment token">#path to your original accurate DEM - make sure path is enclosed by r' '</SPAN>
original_dem <SPAN class="operator token">=</SPAN> r<SPAN class="string token">'c:\\myproject\mygeodatabase\originalDEM'</SPAN>

<SPAN class="comment token">#change the below paths to your interpolated DEM paths </SPAN>
<SPAN class="comment token"># - make sure path is enclosed by r' '</SPAN>
<SPAN class="comment token">#interpolated DEM 1</SPAN>
interp_dem1 <SPAN class="operator token">=</SPAN> r<SPAN class="string token">'c:\\myproject\mygeodatabase\interpDEM1'</SPAN>
<SPAN class="comment token">#interpolated DEM 2</SPAN>
interp_dem2 <SPAN class="operator token">=</SPAN> r<SPAN class="string token">'c:\\myproject\mygeodatabase\interpDEM2'</SPAN>
<SPAN class="comment token">#interpolated DEM 3</SPAN>
interp_dem3 <SPAN class="operator token">=</SPAN> r<SPAN class="string token">'c:\\myproject\mygeodatabase\interpDEM3'</SPAN>
<SPAN class="comment token">#interpolated DEM 4</SPAN>
interp_dem4 <SPAN class="operator token">=</SPAN> r<SPAN class="string token">'c:\\myproject\mygeodatabase\interpDEM4'</SPAN>
<SPAN class="comment token">#interpolated DEM 5</SPAN>
interp_dem5 <SPAN class="operator token">=</SPAN> r<SPAN class="string token">'c:\\myproject\mygeodatabase\interpDEM5'</SPAN>

<SPAN class="comment token">#make a list of the interpolated DEM paths</SPAN>
interp_dem_list <SPAN class="operator token">=</SPAN> <SPAN class="punctuation token">[</SPAN>interp_dem1<SPAN class="punctuation token">,</SPAN> interp_dem2<SPAN class="punctuation token">,</SPAN> interp_dem3<SPAN class="punctuation token">,</SPAN> interp_dem4<SPAN class="punctuation token">,</SPAN> interp_dem5<SPAN class="punctuation token">]</SPAN>

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

<SPAN class="keyword token">print</SPAN><SPAN class="punctuation token">(</SPAN><SPAN class="string token">"Processing Complete"</SPAN><SPAN class="punctuation token">)</SPAN>

<SPAN class="line-numbers-rows"><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN></SPAN>