# Compare two Digital Elevation Model (DEM)

I would like to know what's the best method(s) to compare two DEMs without using any GUI soft like QGIS or ArcGIS.
I'd like to do it in Python or C++, with GDAL or other similar lib.
Does calculating the RMSE (Root mean square error/deviation) is useful?
How to proceed to subtract DEM1 from DEM2? And I should find a flat result, right? Is their other ways?
Because reading for each coordinate the elevation `z` and compare it with the other DEM is a bit heavy I think.

Thanks for help,
eo

• How else would you subtract A from B without comparing each cell? – Nathan W Apr 10 '13 at 8:32
• Don't know, you're right... – eouti Apr 10 '13 at 8:34
• For some methods and insights, you might be interested in reading through a recent case study. Subtracting one DEM from the other (which is quick and easy) is only the very beginning: there will be differences that you have to explore, measure, and seek to understand. Computing the RMSE has its place, but as a single number it's not going to tell you much about how the DEMs differ. – whuber Apr 10 '13 at 15:31
• +1 @whuber. Without knowing much about the DEMs (are they either identical or not identical?) I would think you'd need to consider more than just the difference between matching pixel values. A more in-depth analysis of slope (or other metrics) across a neighborhood would yield more information. To start you can use R for the programmatic calculations and GRASS/GDAL for the file handling without using the GUIs. – Radar Apr 10 '13 at 17:02
• @Radar Yes, slope is important. But you don't need `R` for these calculations: much (maybe most) of what is needed can be accomplished with relatively simple "map algebra" calculations available in ArcGIS/Spatial Analyst or GRASS, for instance. I would recommend `R` when analytical needs are sophisticated and the DEM is relatively small (perhaps a million cells or less), but for large DEMs you need the efficiency of a raster-based GIS (which `R` is definitely not). – whuber Apr 10 '13 at 17:30

``````a = np.random.randint(-10,400,(500,500))      # or gdal.Open("path/to/raster").GetRasterBand(0).ReadAsArray()