I have two land use rasters (6 classes each) and I'd like to create an error matrix (errors of omission and commission) from them both.

The premise is that I have a very fine 'base' raster and have aggregated it to a new, coarser resolution. I would like to know not just the difference in total area of each land use class but from what and to what land class have changed;

# make a raster to simulate fine res base raster;
ras.fine <- raster(nrows=10, ncols=10,xmn=0, xmx=10, ymn=0, ymx=10)
ras.fine[] <- sample(seq(from = 1, to = 6, by = 1), size = 100, replace = TRUE)

# aggregate it to represent the coarser raster
ras.coarse <- aggregate(ras.fine,fact=2,expand=FALSE,fun=modal,na.rm=T)

Once we have done the above, it's easy enough to establish the total areas of land use classes and establish the difference ;

# use freq to count instances of land class and multiply by square of resolution    
data.frame(class = freq(ras.fine)[,1],count = freq(ras.fine)[,2],area = freq(ras.fine)[,2]*(res(ras.fine)[1])^2)

data.frame(class = freq(ras.coarse)[,1],count = freq(ras.coarse)[,2],area = freq(ras.coarse)[,2]*(res(ras.coarse)[1])^2)

but the differences in total area of land class aren't enough; i'd like to know from what land classes the error has occurred from etc, making an error matrix.

As a start, we could directly compare the base raster and the aggregated raster by disaggregating the coarser raster back to the base resolution (is this fair??);

# disaggregate back to original resolution
rc.d <- disaggregate(ras.coarse,fact=2)

# then create a grid of 'disagreement'; that is where the two rasters do not agree
disagree <- ras.fine != rc.d

# and establish what land classes make up those cells of disagreement, from both the fine raster and the coarser raster
fine.cov <- ras.fine * disagree
coarse.cov <- rc.d * disagree

so now I have statistics for how much of each land class is classified incorrectly at a coarser scale, and theoretically what it is going to/coming from.

I'm a bit stuck from here; how do i fashion it into an error matrix? from what land class has another land class commissioned area from and vice versa?

I am essentially analysing the aggregation method but this allows quantification in uncertainty mapped against cost of effort.

Further: Qu on Cross Validation with regards to the validity of error matrices and raster aggregation

2 Answers 2


Why don't you take a sampling approach? Using sampleRandom or sampleRegular with sp = TRUE, you could draw samples from each raster and then just use table. If you used two different sample sizes with sampleRegular you can unalign the sampling grid to revel potential error at representing different scale variation or anisotropy. You could also use spsample with hexagonal sampling, which would be more efficient.

Create example data library(raster) library(sp)

r <- raster(nrows=180, ncols=360, xmn=571823.6, xmx=616763.6, ymn=4423540, 
             ymx=4453690, resolution=300, crs = CRS("+proj=utm +zone=12 +datum=NAD83 
             +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"))
 r[] <- round(runif(ncell(r),1,5),0)
 r2 <- aggregate(r, fact=2, fun=max)

Create extent polygon, create hexagonal sample, extract values and calculate contingency matrix. To extend to area just use the buffer argument in extract and then write a function that calculates area, for each class, that can be passed to lapply. Then you will have two vectors of areas for each sample that can be compared.

ext <- as(extent(r), "SpatialPolygons")
pts <- sp::spsample(ext, type = "hexagonal", cellsize = 900)
  pts <- SpatialPointsDataFrame(pts, data.frame(ID=1:length(pts)))

pts@data <- data.frame(pts@data, r=extract(r, pts))
pts@data <- data.frame(pts@data, r2=extract(r2, pts))

table(pts$r, pts$r2)
  • @ Jeffrey Evans, this is certainly something to consider and explore because my working data is of a very large size and sampling could be a good way to go
    – Sam
    Commented Apr 19, 2016 at 13:55
  • However, i'm not classifying as such, i'm simply assessing the error of aggregation (for many resolutions) from a base resolution. Therefore i know all of the true data.
    – Sam
    Commented Apr 19, 2016 at 13:56
  • Also, another issue i need to think about is clumping, some land classes occur more frequently as bigger areas than others and bigger clumps are less susceptible to error at coarser resolutions than smaller clumps. Therefore could sampling allow for this?
    – Sam
    Commented Apr 19, 2016 at 13:58
  • @Sam, I believe that sampling would be a much more tractable approach. My example is fairly simplistic and you would need to expand it to include the sampling of "area" using the buffer option in extract. This should address your clumping issue. In this context aggregation error and uncertainty are a real issue and the only "true" data is your base resolution raster. If you get stuck on what I am getting at, let me know and I will expand my answer. Although, it may take me a bit to get to it. Commented Apr 19, 2016 at 14:16
  • @ Jeffrey Evans, i'm a bit busy myself but i'll certainly try to apply your methods in the next week or so. i'd like to see what results they produce. Error and uncertainty are, as you say, a real issue and curiously generally overlooked in the rasterization process.
    – Sam
    Commented Apr 19, 2016 at 14:55

actually, i worked something out;

following the line;

rc.d <- disaggregate(ras.coarse,fact=2)

do not create any rasters of disagreement or that sort of thing but stack the base raster and the disaggregated raster and extract the raster cell values in the stack as a data frame, and then shrink this new data frame by counting each unique row;

st <- stack(ras.fine,rc.d)
vdf <- as.data.frame(extract(st,c(1:ncell(st))))
vdf.c <- count(vdf, vars = c("layer.1","layer.2"))

this gives us the 'from' class and the 'to' class and the amount of times it occurs.

Then we can build a master template of all possible changes and merge in our results as extracted from the stack;

template <- data.frame(base=rep(1:6,each=6),new=rep(1:6,6))
results <- merge(template,vdf.c,by.x=c("base","new"),by.y=c("layer.1","layer.2"),all.x=T)
# change NA to zero if required
results$freq[is.na(results$freq)] <- 0

then we can use acast from package reshape2 to quickly turn it into a matrix

mat <- acast(results, base ~ new, value.var='freq')

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.