# Map accuracy assessment by moving window in R

I would like to make an accuracy assessment using a confusion matrix between classified Landsat image and reference dataset. Maps have the same resolution and extent. I would like to evaluate the agreement by pixel-by-pixel. In many studies I have found that they used moving-kernel window (ideal 3x3) to "deal" with the Landsat pixels misregistration. However I couldn't find any approach to do use this moving window for confusion matrix evaluation in R, normally it is used to values interpolation.

Do you have any ideas how to implement moving window into classification accuracy assessment? Or am I misunderstanding this approach?

Thanks a lot,

Example:

library("raster")
clas <- raster(ncols=5, nrows=5)
values(clas)<-c(2,1,0,0,1,
2,1,0,0,0,
1,0,0,0,0,
0,0,0,0,0,
2,1,0,0,2)
reference <- raster(ncols=5, nrows=5)
values(reference)<-c(2,1,0,0,0,
2,1,0,0,1,
0,0,0,0,0,
1,2,0,0,2,
0,0,0,0,0)


EDIT:

After utilisation of movingFun (raster) I think this is mostly to get information in one cell taking into consideration the cells around.

so if I have a

ref<-c(0,0,0,0,1,0,0,0,0)
class<-c(0,0,0,0,1,1,1,0,1)


and I apply movingFun

mov_ref<-movingFun(ref, n=3)
mov_clas<-movingFun(class, n=3)


I will obtain

mov_ref
NA, 0.0000000 0.0000000 0.3333333 0.3333333 0.3333333 0.0000000
0.0000000        NA

mov_clas
NA 0.0000000 0.0000000 0.3333333 0.6666667 1.0000000 0.6666667
0.6666667        NA


For accuracy assessment needs could I directly compare mov_ref and mov_class? Or do I have maybe to reclass it into 0,1 (by threshold 0.5) to compare them? If I leave the numbers as 0.3333 I dont reach a good result in classification accuracy.

• have a look at the documentation of the raster package first. There is for instance a straight-forward function called movingFun() which does exactly what you want. For frequencies of miss-classification you could also look at crosstab – Curlew Jun 29 '14 at 16:45
• Hi, thanks for your response, I have edited my question. The problem is if I can directly use the result from movingFun into Kappa statistic (or crosstab?) thanks a lot – maycca Jun 29 '14 at 18:25

As many on this forum know, I am often for an R solution. However, in this case it is reinventing the wheel, and in a much less robust way. There is a great piece of free software, Map Comparison Kit (MCK), that implements many published and novel validation statistics for rasters. Of particular interest in this case are the Kappa, fuzzy Kappa and weighted Kappa.

Now, if you want to implement something in R there are many approaches you can take that depend on the complexity of the validation statistic. In a univariate case you can easily pass a function to "focal" to calculate uncertainty within a defined neighborhood. Moving into a bivariate case, you would want to vectorize the problem and define a function that would take two independent data into account. I do not believe that "movingFun" or "focal" will take two rasters into account. You can however, use "overlay", "getValuesBlock" or ideally"getValuesFocal" all of which will operate on stack/block objects.

Here is a worked example of calculating Kappa, using a 3x3 window, with "getValuesFocal". In the for loop the lapply function is reclassifying simulated probabilities [p >= t |1| else |0|], The parameter to adjust the sensitivity is "p" and "ws" adjust the size of the focal window extracted. I wrote this to be memory safe so, it writes a file ("Kappa.img") to disk in the defined working directory.

require(raster)
require(asbio)

setwd("D:/TEST")

ws <- 3   # window size
p=0.65    # probability threshold

# Create example data
pred <- raster(ncol=100, nrow=100)
pred[pred] <- runif(length(pred[pred]),0,1)
obs <- pred
obs[obs] <- runif(length(pred[pred]),0,1)
obs.pred <- stack(obs,pred)
names(obs.pred) <- c("obs","pred")

# Create new on-disk raster
s <- writeStart(obs.pred[[1]], "Kappa.img", overwrite=TRUE)
tr <-  blockSize(obs.pred)
options(warn=-1)

# Loop to read raster in blocks using getValuesFocal
for (i in 1:tr$n) { # Get focal values as list matrix object v <- getValuesFocal(obs.pred, row=tr$row[i], nrows=tr$nrows[i], ngb=ws, array=FALSE) # reclassify data to [0,1] using lapply v <- lapply(v, FUN=function(x) { if( length(x[is.na(x)]) == length(x) ) { return( NA ) } else { return( ifelse(x >= p, 1, 0) ) } } ) # Loop to calculate Kappa and assign to new raster using writeValues r <- vector() for( j in 1:dim(v[[1]])[1]) { Obs <- v[[1]][j,] Obs <- Obs[!is.na(Obs)] Pred <- v[[2]][j,] Pred <- Pred[!is.na(Pred)] if( length(Obs) >= 2 && length(Obs) == length(Pred) ) { r <- append(r, Kappa(Pred, Obs)$khat)
} else {
r <- append(r, NA)
}
}
writeValues(s, r, tr\$row[i])
}
s <- writeStop(s)

k <- raster("Kappa.img")
plot(k)


Problem:

• applying moving window in matrix/raster to calculate Kappa statistics between classified map and reference dataset

Solution:

• using focal{raster package}. This will implement moving window in raster
• function "modal" fun=modal to keep the majority values of neighbour values
• movingFun() from the {Raster}is mostly for vectors

R code for focal()

# matrix(characterising moving window) of 3x3pixels
foc.clas<-focal(clas, w=matrix(1,3,3), fun=modal)


In this way I can obtain a new rasters changed by focal function (moving window) from

1. classified map
2. reference dataset

To calculate Kappa statistic, I will use these new-created rasters by focal function. I am using {asbio} package, so firstly the rasters have to be converted as.matrix(clas)

• Could you please add some more detail to your answer? Perhaps a brief worked example. Since you cannot pass a stack to focal, how would you apply a focal function to the "new-created rasters? Recommending matrix coercion of rasters is not the best advice because it will not keep the workflow memory safe. Other than smoothing, I am not clear on how assigning the mode to the focal cell allows a fuzzy assessment of accuracy. You would want to asses all the values in the focal window, which is consistent with the matrix algebra defined in Hagen-Zanker (2006 & 2009). – Jeffrey Evans Nov 13 '14 at 22:16
• Hi Jeffrey, thanks for your answer. However, I dont really understand some parts of your answer: why do you use probability threshold (p=0,65)? - it is only to create binary 0,1 example data? For getValuesFocal I think it is what I needed whan I was asking. My approach firstly converted my rasters: class and ref using focal approach, majority filter to foc.class and foc.ref and after I applied Kappa stat on these rasters. I think your approach is better. Thanks also for MCK tool. – maycca Dec 5 '14 at 9:37
• I used a p threshold to make the answer a bit broader in application. Often a binary response is derived from an estimated probability distribution. This is only relevant to the block of code where I am creating example data and does not relate to the implementation of the moving window Kappa. If my answer works for you please marked answered so we can close this thread. – Jeffrey Evans Dec 5 '14 at 16:26