# Comparing two Digital Elevation Models (DEMs) from LAS files?

I have two Lidar las files, one is original let's say with X points. And the other is copy of the first las file but with Y points, where Y is less than X.

Now, I want to compare how the Digital Elevation Models of these two las files vary...

I want to get information like RMSE, standard deviation, etc...

I would appreciate, if anyone could tell me what softwares, or ways to get the comparison info...

-
Perhaps you could edit the question to use M and N (instead of X and Y). At first reading I thought X and Y were the coordinate values! – Mark Ireland Oct 14 '10 at 17:15
You really do need to provide more information to receive relevant help. Your current question makes very little scene. Which field in your las file is holding the values. The way ground classified are assigned in the las format are a classification field and not different z (elevation) values. A vendor would have had to used on of the unassigned fields to hold a difference in z values. – Jeffrey Evans Apr 7 '13 at 23:21

How to compare two Digital Elevation Models (DEMs)

Solution using the software R.

#-------------------------------------------------------------------------
#Creating a reproducible example

library(raster)

#simulating raster_1

f = system.file("external/test.grd", package="raster")
DEM_1 = raster(f)

#simulating raster_2

DEM_2 = DEM_1
# replacing values from raster_1 to create a new raster sample (raster_2)
DEM_2[(DEM_2>500 & DEM_2<900)] = 550
DEM_2[(DEM_2>200 & DEM_2<300)] = 500

#-------------------------------------------------------------------------
# Comparison 1 (DEM_3 resulted from subtracting DEM_2 from DEM_1)

DEM_3 = DEM_1 - DEM_2

par(mfrow=c(1,3))

plot(DEM_1, main = "DEM_1")
plot(DEM_2, main = "DEM_2")
plot(DEM_3, main = "DEM_3 = DEM_1 - DEM_2")

dev.off()

#-------------------------------------------------------------------------
#Comparison 2 (histogram)

hist(DEM_1, prob=T, main="DEM_1", xlab="")
hist(DEM_2, prob=T, main="DEM_2", xlab="")
hist(DEM_3, prob=T, main="DEM_3 = DEM_1 - DEM_2", xlab="")

par(mfrow=c(1,1))

standard_deviation = sd(c(as.matrix(DEM_3)),na.rm=T)

dev.off()

#-------------------------------------------------------------------------
#comparison 3 (RMSE)

library(hydroGOF)

DEM_1_matrix = c(as.matrix(DEM_1))
DEM_2_matrix = c(as.matrix(DEM_2))

rmse = rmse(DEM_1_matrix,DEM_2_matrix)
rmse
[1] 135.3675 # this is the root mean squared error (RMSE) result.

The second and third parts of previous version of this answer were improved and migrated to a more related question, here. This edit was made aiming to improve readability on both answers.

-

You can do this through the ESRI ArcGIS Geostatistical Analysis Extension - there is a section in the help on performing validation on subsets.

You could do the same through GRASS through the R interface. Tomislav Hengl describes in some detail how to do so in his book A Practical Guide to Geostatistical Mapping. It's open access, so the PDF is free to download.

-

I would do a simple DEM of difference. DEM2-DEM1. This will show all areas that are different and by how much.

Theres an image to a high res dem of difference on my website homepage. thadwester.com
Take a look at the colorful left image.

-

As far as I know, RMSE is only stated during the making of the DEM, and not as an attribute for further refrence, so you'de have to "catch it" manually during the making of the DEM (that said, I never made a DEM from Lidar, only from other data).

If you want to see the differences between the DATA inside the two DEMS, I'd use cut/fill which is in the Spatial Analyst extension of ArcGIS (under "Surface Analysis"). The cut/fill shows you in a simple thematic map the changes between the DEM's.

-
Cut/fill is too crude because it does not quantify the differences. RMSE is a general way to compare two datasets: not only is it useful to compare a DEM to ground-truth data, it is one way to quantify differences between two DEMs. – whuber Oct 14 '10 at 14:54