Graphically depict statistical difference between thematic rasters

Question

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)

# Read the TIF data
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[1], 2)
b1 <- round(coefs[2],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)[2], range(old)[2])
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

Answer 1

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!