# Graphically depict statistical difference between thematic rasters

The graph shows a regression between two raster datasets with classified tree canopy cover: class 1 = 0-10%, class 2 = 10 - 20% ...class 6 = >60% (Figure 1, Figure 2). The spatial resolution for raster 1 is 1m while the spatial resolution for raster 2 is 30m. Note that raster 2 was resampled to 1m for analysis purposes in the script (Appendix A).

I am looking for an elegant, graphical approach to statistically describe the difference between two thematic rasters at different spatial scales. My first attempt was to 1) create a raster stack, 2) randomly sample the stack, 3) run a regression and descriptive statistics, 4) graph the results (Appendix A). What other graphical methods exist for statistically comparing thematic raster data? I usually use R, Python and MATLAB, although I would be happy with generic solutions too.

Figure 1 Figure 2 Appendix A

``````require(raster)
require(spatstat)
require(ggplot2)
require(gridExtra)
require(hexbin)
require(rgdal)

file = 'C:/path/raster1.tif'
file2 = 'C:/path/raster2.tif'

# Create raster objects from the single band canopy products
r = raster(file)
r2 = raster(file2)

# Resample r2 to 1m for analysis purposes
r2 = resample(r2, r, method = "ngb")

# Stack the two raster layers for analysis
raster <- stack(r, r2)

# sample random locations in raster stack and report values in a dataframe
df = data.frame(sampleRandom(raster, size=1000, cells=TRUE, sp=TRUE))

## Do the regression and plot the results
new = df\$X3711203_ne
old = df\$X3711203_ne_30m

# Calculate RMSE and other values
fit <- lm(old ~ new)
rmse <- round(sqrt(mean(resid(fit)^2)), 2)
coefs <- coef(fit)
b0 <- round(coefs, 2)
b1 <- round(coefs,2)
r <- round(sqrt(summary(fit)\$r.squared), 2)

# Build equation see ?plotmath
eqn <- bquote(italic(y) == .(b0) + .(b1)*
italic(x) * "," ~~ r == .(r) * "," ~~ RMSE == .(rmse))

eqn_text = as.character(as.expression(eqn))

# Get max range to automatically update plot size
li = c(range(new), range(old))
range = tail(sort(li),1) + 0.05

# Plot results
p1 = ggplot(df, aes(new, old))  +
geom_point() +
geom_smooth(method = lm) +
stat_sum(aes(size = ..n..)) +
xlab("Cover Class (1m)") +
ylab("Cover Class (30m)") +
annotate("text", x= 0, y= range, label=eqn_text, hjust=0, size=4, face="italic", parse=TRUE) +
xlim(0, range) +
ylim(0, range) +
ggtitle("Tree Canopy Cover Class")

p1
``````
• Generally linear regression is used when when the predictor and response variables are continuous. Your variables are categorical. You could try an ANOVA analysis, by comparing the mean area of each cover category at your different resolutions. That would make your predictor variable (cover%) categorical and your response (area) continuous. – khafen Sep 3 '14 at 20:51

You can also make a raster of "agreement".

One raster can be binary 0=disagreement and 1=agreement.

And another raster depicting how far off you are...for example, if a certain pixel has been classified as 1 but is actually a 3, then your difference is 2. In your case, the difference scale would go from 0 to 6...colour-code it differently than your actual rasters to avoid confusion and you end up with a nice visual showing degree of agreement with spatial context.

With python+numpy, it would be pretty easy.

Binary raster would be:

``````binary_agreement_raster = (predicted_raster == truth_raster).astype(np.int)
``````

and the degree of agreement:

``````degree_agreement = np.abs(predicted_raster - truth_raster)
``````

good to go!