How to calculate RMSE, skewness, and kurtosis of raster image?

Question

I'm doing a project to analyze the impact different interpolation techniques plays on DEM accuracy.

In order to do so, I need to calculate the RMSE, skewness, and kurtosis.

I have:

1. an original source DEM,
2. a DEM generated from 10 different interpolation techniques,
3. a DEM generated by the subtraction of above two DEMs.

How do I calculate the three parameters?

Answer 1

If you have a statistically sufficient number of survey data points you could use the elevation values in those points to calculate the RSME using the code below. The survey data table would need to have some existing fields including the surveyed elevation, and the interpolated elevations (from the add surface information tool in the 3d analyst toolbox).

Here is the code written for 9.3.

try:

import arcgisscripting, exceptions, math, traceback, sys
gp = arcgisscripting.create(9.3)
gp.OverWriteOutput = 1
gp.toolbox = "management"



TheFeatureClassOfPoints = "Z:\Data\Cadastral\cadastrepoints\Cadastrepoints.gdb\Points\LummiDeltaElevationSurveyControlPoints"
InterpolatedValueField = r"deltaz"
SurveyedValue = r"surveyz"
RMSEField = r"RMSEdelta"


SquaredValues = []


print "Start cursor"
rows = gp.SearchCursor(TheFeatureClassOfPoints)
row = rows.next()
while row:
    i = row.GetValue(InterpolatedValueField)
    s = row.GetValue(SurveyedValue)
    SquaredValues.append(pow((i - s), 2))
    row = rows.next()

del row, rows    

n = len(SquaredValues)

SumOfSquares = sum(SquaredValues)

x = SumOfSquares/n

rmse = str(pow(x, .5))
print "rmse =", rmse
print "Add field."

gp.AddField_management(TheFeatureClassOfPoints, RMSEField, "DOUBLE")
print "Calculating field."
gp.CalculateField(TheFeatureClassOfPoints, RMSEField, rmse, "VB")

del gp
print "Finished"

Geoprocessing Errors will be caught here

except arcgisscripting.ExecuteError: print gp.GetMessages(2) raw_input("ArcGIS Error, press enter to continue")

Other errors will get caught here

except Exception, ErrorDesc: print ErrorDesc.message import traceback, sys # get the traceback object tb = sys.exc_info()[2]

# tbinfo contains the errors line number and the code   
tbinfo = traceback.format_tb(tb)[0]

print tbinfo        # provides where the error occurred
print sys.exc_type  # provides the type of error
print sys.exc_value # provides the error message
raw_input("Python Error, press enter to continue")

Answer 2

You can do this in R. Here's a script that reads in your DEMs, calculates your stats, plots your results and then writes to csv (if you really want to plot in a spreadsheet). Remember to add your own file locations and you may need to tweak the RMSE, I didn't check that.

install.packages("raster")
install.packages("e1071")
library(raster)
library(e1071)

# your original raster
org = raster("~/dir/dir/file.tif")
org = na.omit(getValues(org))

# your computed rasters
f = list.files("~/dir/dir", full.names = T, pattern = "tif")
# for naming results
n = list.files("~/dir/dir", pattern = "tif")
n = lapply(f, function(i){unlist(strsplit(i, "[.]"))[1]})
# read computed rasters
r = lapply(f, raster)
names(r) = n
# convert to a numeric vector
v = lapply(r, function(i){na.omit(getValues(i))})

# Calculate your values
kurt = sapply(v, kurtosis)
skew = sapply(v, skewness)
# RMSE
RMSE = sapply(v, function(i){
   sqrt(mean((i - org)^2, na.rm = TRUE))
})

# Wrap up your results
df = data.frame(DEM = names(v), kurt, skew, RMSE)

# Plot results
barplot(kurt)
# etc

# Write results
write.csv(df, "~/dir/dir/file.csv", rownames=F)

Answer 3

I calculated the rmse by importing the excel sheet in matlab...thank you all,I'm still on the process by adapting new procedures.

Answer 4

You can do this using GRASS, providing your DEMs are integer (a restrictive limitation). You need the r.statistics tool and also r.mapcalc (for RMSE).