3

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

Improve this question
11
  • How else would you subtract A from B without comparing each cell?
    – Nathan W
    Commented Apr 10, 2013 at 8:32
  • 1
    Don't know, you're right...
    – eouti
    Commented Apr 10, 2013 at 8:34
  • 2
    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
    Commented Apr 10, 2013 at 15:31
  • 1
    +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
    Commented Apr 10, 2013 at 17:02
  • 1
    @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
    Commented Apr 10, 2013 at 17:30

1 Answer 1

Reset to default
3

There's a similar question, but not doing it with a programming language.

here's a small python snippet which checks the difference between two numpy matrixes.

a = np.random.randint(-10,400,(500,500))      # or gdal.Open("path/to/raster").GetRasterBand(0).ReadAsArray()
b = np.random.randint(-10,400,(500,500))
dif = abs((a-b)/(b+1e-5))<0.05                # Returns a T/F matrix if difference is more than 5%. 
                                              # I add 1e-5 at the denominator to avoid div by zero
float(np.count_nonzero(tf))/len(tf.flatten()) # Ratio between Trues vs length of matrix
Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.