In colinr23/spatialcompare: Spatial Pattern Comparison Metrics

knitr::opts_chunk$set(echo = TRUE)

Area Overlap Statistic

The area overlap statistics include all the statistics described by @maruca_area-based_2002. These statistics are probabilistic pattern association tests that are appropriate when edge effects are present, polygon size is heterogeneous, and the number of polygons varies from one classification to another.

Relative Area Overlap

The relative area overlap is area of intersection as a fraction of the area of union of two polygons. For identical polygons the relative area overlap will be 1.

library(gpclib)
library(sp)
library(spatialcompare)

set.seed(1234)
theta <- seq(0, 2 * pi, length=(100))
poly1<-  cbind(c(0 + 1 * cos(theta) + rnorm(100, sd=0.1)), c(0 + 2 * sin(theta)))
poly2<-  cbind(c(0 + 2 * cos(theta) ), c(-1 + 1.5 * sin(theta)+ rnorm(100, sd=0.1)))
x = 300
y = 300
plot(x, y, type = "n", xlim=c(-5,5), ylim=c(-5,5))
polygon(poly1)
polygon(poly2)
p1 <- as(poly1, "gpc.poly")
p2 <- as(poly2, "gpc.poly")
relative_area_overlap(p1,p2)

Maximum Relative Area Overlap

The maximum relative area overlap gives the the relative area overlap for the polygon which has maximum overlap with the candidate polygon. This can be useful when for example, multiple polygons overlap at the edges but a central polygon overlaps a much larger portion. This is common in image classification and segmentation problems and change analysis.

square <- t(replicate(50, {
  o <- runif(2)
  c(o, o + c(0, 0.1), o + 0.1, o + c(0.1, 0), o)
}))
square2 <- t(replicate(50, {
  o <- runif(2)
  c(o, o + c(0, 0.2), o + 0.2, o + c(0.2, 0), o)
}))
ID <- paste0('sq', seq_len(nrow(square)))
ID2 <- paste0('sq', seq_len(nrow(square2)))
#Create SP
polys <- SpatialPolygons(mapply(function(poly, id) {
  xy <- matrix(poly, ncol=2, byrow=TRUE)
  Polygons(list(Polygon(xy)), ID=id)
}, split(square, row(square)), ID))
polys2 <- SpatialPolygons(mapply(function(poly, id) {
  xy <- matrix(poly, ncol=2, byrow=TRUE)
  Polygons(list(Polygon(xy)), ID=id)
}, split(square2, row(square2)), ID2))
# Create SPDF
polys.df <- SpatialPolygonsDataFrame(polys, data.frame(id=ID, row.names=ID))
polys.df2 <- SpatialPolygonsDataFrame(polys2, data.frame(id=ID, row.names=ID))
plot(polys.df, col=rainbow(50, alpha=0.5))
ids <- sapply(polys.df@polygons, function(p) p@ID)
plot(polys.df2, col=rainbow(50, alpha=0.5))
ids2 <- sapply(polys.df2@polygons, function(p) p@ID)

p <- polys2@polygons[[1]]@Polygons[[1]]@coords  #Sample Polygon for maximum relative area Overlap

maximum_rel_area_ovelap(polys,p)

Average Maximum Area Overlap

This function averages the maximum area overlap over the entire dataset. This is called the area overlap statistic in [maruca_area-based_2002].

average_max_rel_area_overlap(polys,polys2)

Simultaneous Area Overlap

Since overlap can be evaluated two ways (I:J and J:I) and is asymmetric, the simultaneous version gives a general overlap statistic for the two sets of polygons.

simultaneous_area_overlap(polys,polys2)

Weighted Average Relative Area Overlap

The average maximum relative overlap can also be computed as a weighted average where weights are given by the area of the focus polygon.

weighted_avg_maximum_relative_area_overlap(polys,polys2)

Weighted Simultaneous Area Overlap

weighted_simultanous_area_overlap(polys,polys2)

Sample Analysis

Data from [@maruca_area-based_2002] were re-created to demonstrate metric values for overlap statistics.

  data('scenarios')
  proj4string(scenarios) <- ""
  Scenes = vector(length=10)
  Ai = vector(length=10)
  Aj = vector(length=10)
  Aij = vector(length=10)

  #scenario 1 - perfect overlap
  Scenes[1] <- "1. Perfect overlap-unweighted"
  Scenes[2] <- "1. Perfect overlap-weighted"
  #area overlap - I
  Ai[1] <- average_max_rel_area_overlap(subset(scenarios, Scenario == 1 & Poly == 1), subset(scenarios, Scenario == 1 & Poly == 1))
  Ai[2] <- weighted_avg_maximum_relative_area_overlap(subset(scenarios, Scenario == 1 & Poly == 1), subset(scenarios, Scenario == 1 & Poly == 1))
  #area overlap - J
  Aj[1] <- average_max_rel_area_overlap(subset(scenarios, Scenario == 1 & Poly == 1), subset(scenarios, Scenario == 1 & Poly == 1))
  Aj[2] <- weighted_avg_maximum_relative_area_overlap(subset(scenarios, Scenario == 1 & Poly == 1), subset(scenarios, Scenario == 1 & Poly == 1))
  #area overlap - simultanious
  Aij[1] <- simultaneous_area_overlap(subset(scenarios, Scenario == 1 & Poly == 1), subset(scenarios, Scenario == 1 & Poly == 1))
  Aij[2] <- weighted_simultanous_area_overlap(subset(scenarios, Scenario == 1 & Poly == 1), subset(scenarios, Scenario == 1 & Poly == 1))

  #scenario 2 - 
 #area overlap - I
  Scenes[3] <- "2. Maximally offset-unweighted"
  Scenes[4] <- "2. Maximally offset-weighted"
  Ai[3] <- average_max_rel_area_overlap(subset(scenarios, Scenario == 2 & Poly == 1), subset(scenarios, Scenario == 2 & Poly == 2))
  Ai[4] <- weighted_avg_maximum_relative_area_overlap(subset(scenarios, Scenario == 2 & Poly == 1), subset(scenarios, Scenario == 2 & Poly == 2))
  #area overlap - J
  Aj[3] <- average_max_rel_area_overlap(subset(scenarios, Scenario == 2 & Poly == 2), subset(scenarios, Scenario == 2 & Poly == 1))
  Aj[4] <- weighted_avg_maximum_relative_area_overlap(subset(scenarios, Scenario == 2 & Poly == 2), subset(scenarios, Scenario == 2 & Poly == 1))
  #area overlap - simultanious
  Aij[3] <- simultaneous_area_overlap(subset(scenarios, Scenario == 2 & Poly == 1), subset(scenarios, Scenario == 2 & Poly == 2))
  Aij[4] <- weighted_simultanous_area_overlap(subset(scenarios, Scenario == 2 & Poly == 1), subset(scenarios, Scenario == 2 & Poly == 2))

    #scenario 3 - 
 #area overlap - I
  Scenes[5] <- "3. Good overlap - polygon areas similar -unweighted"
  Scenes[6] <- "3. Good overlap - polygon areas similar -weighted"
  Ai[5] <- average_max_rel_area_overlap(subset(scenarios, Scenario == 3 & Poly == 1), subset(scenarios, Scenario == 3 & Poly == 2))
  Ai[6] <- weighted_avg_maximum_relative_area_overlap(subset(scenarios, Scenario == 3 & Poly == 1), subset(scenarios, Scenario == 3 & Poly == 2))
  #area overlap - J
  Aj[5] <- average_max_rel_area_overlap(subset(scenarios, Scenario == 3 & Poly == 2), subset(scenarios, Scenario == 3 & Poly == 1))
  Aj[6] <- weighted_avg_maximum_relative_area_overlap(subset(scenarios, Scenario == 3 & Poly == 2), subset(scenarios, Scenario == 3 & Poly == 1))
  #area overlap - simultanious
  Aij[5] <- simultaneous_area_overlap(subset(scenarios, Scenario == 3 & Poly == 1), subset(scenarios, Scenario == 3 & Poly == 2))
  Aij[6] <- weighted_simultanous_area_overlap(subset(scenarios, Scenario == 3 & Poly == 1), subset(scenarios, Scenario == 3 & Poly == 2))

    #scenario 4 - 
  Scenes[7] <- "4. Good overlap - polygon areas different -unweighted"
  Scenes[8] <- "4. Good overlap - polygon areas different -weighted"
 #area overlap - I
  Ai[7] <- average_max_rel_area_overlap(subset(scenarios, Scenario == 4 & Poly == 1), subset(scenarios, Scenario == 4 & Poly == 2))
  Ai[8] <- weighted_avg_maximum_relative_area_overlap(subset(scenarios, Scenario == 4 & Poly == 1), subset(scenarios, Scenario == 4 & Poly == 2))
  #area overlap - J
  Aj[7] <- average_max_rel_area_overlap(subset(scenarios, Scenario == 4 & Poly == 2), subset(scenarios, Scenario == 4 & Poly == 1))
  Aj[8] <- weighted_avg_maximum_relative_area_overlap(subset(scenarios, Scenario == 4 & Poly == 2), subset(scenarios, Scenario == 4 & Poly == 1))
  #area overlap - simultanious
  Aij[7] <- simultaneous_area_overlap(subset(scenarios, Scenario == 4 & Poly == 1), subset(scenarios, Scenario == 4 & Poly == 2))
  Aij[8] <- weighted_simultanous_area_overlap(subset(scenarios, Scenario == 4 & Poly == 1), subset(scenarios, Scenario == 4 & Poly == 2))

    #scenario 5 - 
  Scenes[9] <- "5. Overlapping partitions - different spatial scales -unweighted"
  Scenes[10] <- "5. Overlapping partitions - different spatial scales -weighted"
 #area overlap - I
  Ai[9] <- average_max_rel_area_overlap(subset(scenarios, Scenario == 5 & Poly == 1), subset(scenarios, Scenario == 5 & Poly == 2))
  Ai[10] <- weighted_avg_maximum_relative_area_overlap(subset(scenarios, Scenario == 5 & Poly == 1), subset(scenarios, Scenario == 5 & Poly == 2))
  #area overlap - J
  Aj[9] <- average_max_rel_area_overlap(subset(scenarios, Scenario == 5 & Poly == 2), subset(scenarios, Scenario == 5 & Poly == 1))
  Aj[10] <- weighted_avg_maximum_relative_area_overlap(subset(scenarios, Scenario == 5 & Poly == 2), subset(scenarios, Scenario == 5 & Poly == 1))
  #area overlap - simultanious
  Aij[9] <- simultaneous_area_overlap(subset(scenarios, Scenario == 5 & Poly == 1), subset(scenarios, Scenario == 5 & Poly == 2))
  Aij[10] <- weighted_simultanous_area_overlap(subset(scenarios, Scenario == 5 & Poly == 1), subset(scenarios, Scenario == 5 & Poly == 2))
#NOTE i and j appeared to be mixed up. Reversing corrects it.. 
df <- data.frame(Scenario = Scenes, Ai_Value = Aj, Aj_Value = Ai, Aij_Value = Aij)

Here is some fake data replicating the four offset scenarios (scenarios 2-5) in [@maruca_area-based_2002] and a table giving the values of the statistics computing using the spatialcompare package.

par(mar = c(0.001, 0.1, 0.8, 0.1), mfrow=c(2,2))
plot(subset(scenarios, Scenario == 2 & Poly == 1), col="grey", main="Scenario 2")
plot(subset(scenarios, Scenario == 2 & Poly == 2), add=TRUE, border="green")
plot(subset(scenarios, Scenario == 3 & Poly == 1), col="grey", main="Scenario 3")
plot(subset(scenarios, Scenario == 3 & Poly == 2), add=TRUE, border="green")
plot(subset(scenarios, Scenario == 4 & Poly == 1), col="grey", main="Scenario 4")
plot(subset(scenarios, Scenario == 4 & Poly == 2), add=TRUE, border="green")
plot(subset(scenarios, Scenario == 5 & Poly == 1), col="grey", main="Scenario 5")
plot(subset(scenarios, Scenario == 5 & Poly == 2), add=TRUE, border="green")

knitr::kable(df)

Quantity Disagreement

The Quantity Disgreement is a statistic introduced in [jr_death_2011] ("http://www.tandfonline.com/doi/abs/10.1080/01431161.2011.552923"). This statistic along with the allocation disagreement statistic are presented as an alternative to the kappa indices in the paper.

Setting up the data

A1 = matrix(c(0, 0, 0, 0, 0, 0,0,0,0),nrow=3, ncol=3, byrow = TRUE)
A2 = matrix(c(0, 0, 0, 0, 0, 0,0,0,0),nrow=3, ncol=3, byrow = TRUE)
refer = matrix(c(0, 0, 0, 0, 0, 0,1,1,1),nrow=3, ncol=3, byrow = TRUE)
comparison = matrix(c(0, 0, 1, 0, 0, 0,1,0,0),nrow=3, ncol=3, byrow = TRUE)
A3 = matrix(c(0,1,1,0,0,0,1,0,0), nrow = 3, ncol = 3, byrow = TRUE)
A4 = matrix(c(0,0,0,1,0,0,1,1,1),nrow = 3, ncol = 3, byrow = TRUE)

Quantity Disagreement

quantity_disagreement(comparison,refer)
quantity_disagreement(A1,refer)
quantity_disagreement(A2,refer)
quantity_disagreement(A3,refer)
quantity_disagreement(A4,refer)

Allocation Disagreement

The Allocation Disgreement is a statistic introduced in the paper Death to Kappa: birth of quantity disagreement and allocation disagreement for accuracy assessment. This statistic along with the quantity disagreement statistic are presented as an alternative to the kappa indices in the paper.

allocation_disagreement(comparison,refer)
allocation_disagreement(A1,refer)
allocation_disagreement(A2,refer)
allocation_disagreement(A3,refer)
allocation_disagreement(A4,refer)

Structural Similarity Index

The SSIM is a statistic introduced in the paper "Image Quality Assessment: From Error Visibility to Structural Similarity" available here. This is a method used to assessing perceptual image quality.

data("einsteins")
library(raster)
par(mfrow=c(2,3))
image(einsteins[[1]], col=grey(seq(0, 1, length = 256)), axes=FALSE, main=paste("sssim =", round(ssim(einsteins[[1]],einsteins[[1]],w = 5)[[1]], 3)))
image(einsteins[[2]], col=grey(seq(0, 1, length = 256)), axes=FALSE, main=paste("sssim =", round(ssim(einsteins[[1]],einsteins[[2]],w = 5)[[1]], 3)))
image(einsteins[[3]], col=grey(seq(0, 1, length = 256)), axes=FALSE, main=paste("sssim =", round(ssim(einsteins[[1]],einsteins[[3]],w = 5)[[1]], 3)))
image(einsteins[[4]], col=grey(seq(0, 1, length = 256)), axes=FALSE, main=paste("sssim =", round(ssim(einsteins[[1]],einsteins[[4]],w = 5)[[1]], 3)))
image(einsteins[[5]], col=grey(seq(0, 1, length = 256)), axes=FALSE, main=paste("sssim =", round(ssim(einsteins[[1]],einsteins[[5]],w = 5)[[1]], 3)))
image(einsteins[[6]], col=grey(seq(0, 1, length = 256)), axes=FALSE, main=paste("sssim =", round(ssim(einsteins[[1]],einsteins[[6]],w = 5)[[1]], 3)))