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.
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, 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