Skip to content

Commit

Permalink
Include xtable method for twdtw assessment
Browse files Browse the repository at this point in the history
  • Loading branch information
vwmaus committed Feb 7, 2017
1 parent 17e4eec commit 59ece52
Show file tree
Hide file tree
Showing 43 changed files with 2,595 additions and 23 deletions.
5 changes: 3 additions & 2 deletions DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,8 @@ Imports:
sp,
lubridate,
caret,
mgcv
mgcv,
xtable
Suggests:
testthat,
knitr,
Expand All @@ -45,7 +46,6 @@ Suggests:
grid,
png,
Hmisc,
xtable,
tikzDevice
License: GPL (>= 2) | file LICENSE
URL: https://github.com/vwmaus/dtwSat/
Expand Down Expand Up @@ -90,5 +90,6 @@ Collate:
'twdtwAssess.R'
'twdtwClassify.R'
'twdtwCrossValidation.R'
'xtable.R'
'zzz.R'
VignetteBuilder: knitr
3 changes: 3 additions & 0 deletions NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -67,6 +67,7 @@ exportMethods(twdtwCrossValidation)
exportMethods(twdtwMatches)
exportMethods(twdtwRaster)
exportMethods(twdtwTimeSeries)
exportMethods(xtable)
import(ggplot2)
import(methods)
import(raster)
Expand Down Expand Up @@ -110,6 +111,8 @@ importFrom(stats,qnorm)
importFrom(stats,sd)
importFrom(stats,window)
importFrom(stats,xtabs)
importFrom(xtable,print.xtable)
importFrom(xtable,xtable)
useDynLib(dtwSat,bestmatches)
useDynLib(dtwSat,computecost)
useDynLib(dtwSat,g)
Expand Down
6 changes: 4 additions & 2 deletions R/methods.R
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,6 @@
# #
###############################################################


setGeneric("layers",
function(x) standardGeneric("layers"))

Expand Down Expand Up @@ -455,6 +454,8 @@ show.twdtwAssessment = function(object){
print(object@accuracySummary$UsersAccuracy, digits=2)
cat("\nProducers\n")
print(object@accuracySummary$ProducersAccuracy, digits=2)
cat("\nArea and uncertainty\n")
print(object@accuracySummary$AreaUncertainty, digits=2)
invisible(NULL)
}

Expand Down Expand Up @@ -495,7 +496,7 @@ summary.twdtwCrossValidation = function(object, conf.int=.95, ...){
names(l_names) = l_names
ic_ov = mean_cl_boot(x = ov[, c("OV")], conf.int = conf.int, ...)
names(ic_ov) = NULL
assess_ov = data.frame(Accuracy=ic_ov[1], sd=sd_ov, CImin=ic_ov[2], CImax=ic_ov[3])
assess_ov = unlist(c(Accuracy=ic_ov[1], sd=sd_ov, CImin=ic_ov[2], CImax=ic_ov[3]))
ic_ua = t(sapply(l_names, function(i) mean_cl_boot(x = uapa$UA[uapa$label==i], conf.int = conf.int, ...)))
names(ic_ua) = NULL
assess_ua = data.frame(Accuracy=unlist(ic_ua[,1]), sd=sd_uapa[,"UA"], CImin=unlist(ic_ua[,2]), CImax=unlist(ic_ua[,3]))
Expand Down Expand Up @@ -577,3 +578,4 @@ setMethod("projecttwdtwRaster", "twdtwRaster",




41 changes: 23 additions & 18 deletions R/twdtwAssess.R
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ setGeneric("twdtwAssess",
#' @inheritParams twdtwAssessment-class
#' @aliases twdtwAssess
#'
#' @describeIn This function performs an accuracy assessment
#' @describeIn twdtwAssessment This function performs an accuracy assessment
#' of the classified maps. The function returns Overall Accuracy,
#' User's Accuracy, Produce's Accuracy, error matrix (confusion matrix),
#' and estimated area according to [1]. The function returns the metrics
Expand Down Expand Up @@ -54,7 +54,12 @@ setGeneric("twdtwAssess",
#' # Assess classification
#' twdtw_assess = twdtwAssess(r_lucc, validation_samples, proj4string=proj_str)
#' twdtw_assess
#'
#'
#'
#' xtable(twdtw_assess, type="matrix")
#' xtable(twdtw_assess, type="accuracy")
#' xtable(twdtw_assess, type="area")
#'
#' }
#' @export
setMethod(f = "twdtwAssess", signature = "twdtwRaster",
Expand Down Expand Up @@ -134,7 +139,7 @@ twdtwAssess.twdtwRaster = function(object, y, labels, id.labels, proj4string, co

# Error matrix

error_matrix = data.frame(cbind(x, Total=total_map, A=mapped_area, w=w))
error_matrix = data.frame(cbind(x, Total=total_map, Area=mapped_area, w=w))
error_matrix["Total",] = colSums(error_matrix)

# Proportions
Expand All @@ -144,7 +149,7 @@ twdtwAssess.twdtwRaster = function(object, y, labels, id.labels, proj4string, co
total_prop_ref = colSums(y, na.rm = TRUE)

# Proportions matrix
prop_matrix = data.frame(y, Total = total_prop_map, A = mapped_area, w = w)
prop_matrix = data.frame(y, Total = total_prop_map, Area = mapped_area, w = w)
prop_matrix["Total",] = colSums(prop_matrix, na.rm = TRUE)

# Accuracy
Expand All @@ -157,9 +162,6 @@ twdtwAssess.twdtwRaster = function(object, y, labels, id.labels, proj4string, co
names(PA) = cnames
OA = sum(diag(as.matrix(prop_matrix[cnames,cnames])), na.rm = TRUE)

#a_pixel = as.numeric(prop_matrix["Total",cnames] * prop_matrix["Total","A"])
#names(a_pixel) = cnames
#a_ha = a_pixel*res^2/100^2
temp = w^2*UA*(1-UA)/(total_map-1)

VO = sum(temp, na.rm = TRUE)
Expand All @@ -170,38 +172,41 @@ twdtwAssess.twdtwRaster = function(object, y, labels, id.labels, proj4string, co
SU = sqrt(VU)
UCI = SU * mult

fun1 = function(x, xt, A){
sum(A*x/xt, na.rm = TRUE)
fun1 = function(x, xt, Area){
sum(Area*x/xt, na.rm = TRUE)
}

fun2 = function(i, x, xt, A, PA){
fun2 = function(i, x, xt, Area, PA){
x = as.numeric(x[,i])
x = x[-i]
xt = xt[-i]
A = A[-i]
Area = Area[-i]
PA = PA[i]
PA^2*sum(A^2*x/xt*(1-x/xt)/(xt-1), na.rm = TRUE)
PA^2*sum(Area^2*x/xt*(1-x/xt)/(xt-1), na.rm = TRUE)
}

Nj = apply(x, 2, fun1, total_map, mapped_area)
expr1 = mapped_area^2*(1-PA)^2*UA*(1-UA)/(total_map-1)
expr2 = sapply(1:nrow(x), fun2, x=x, xt=total_map, A=mapped_area, PA=PA)
# VP = (1/Nj^2)*(expr1+expr2)
expr2 = sapply(1:nrow(x), fun2, x=x, xt=total_map, Area=mapped_area, PA=PA)
VP = (1/sapply(Nj, function(x) ifelse(x==0, 1, x))^2)*(expr1+expr2)
SP = sapply(VP, function(x) ifelse(x==0, 0, sqrt(x)))
PCI = SP * mult

# Compute adjusted area
estimated_area = prop_matrix["Total",cnames] * prop_matrix["Total","A"]
estimated_area = prop_matrix["Total",cnames] * prop_matrix["Total","Area"]
sd_error = apply(prop_matrix[cnames,cnames], 2, function(x) sqrt(sum( (prop_matrix[cnames,"w"]*x[cnames]-x[cnames]^2)/(error_matrix[cnames,"Total"]-1) )) )
sd_error_estimated_area = sd_error * prop_matrix["Total","A"]
sd_error_estimated_area = sd_error * prop_matrix["Total","Area"]
CI_estimated_area = sd_error_estimated_area * mult

res = list(OverallAccuracy = c(Accuracy=OA, Var=VO, sd=SO, ci=OCI),
UsersAccuracy = cbind(Accuracy=UA, Var=VU, sd=SU, ci=UCI),
ProducersAccuracy = cbind(Accuracy=PA, Var=VP, sd=SP, ci=PCI),
EstimateArea = cbind(Mapped=c(prop_matrix[cnames,"A"]), Estimated=c(estimated_area), ci=c(CI_estimated_area)),
ErrorMatrix = error_matrix
AreaUncertainty = cbind(Mapped=c(prop_matrix[cnames,"Area"]),
Adjusted=c(estimated_area),
ci=c(CI_estimated_area)),
ErrorMatrix = error_matrix,
ProportionMatrix = prop_matrix,
conf.int = conf.int
)

res
Expand Down
2 changes: 1 addition & 1 deletion R/twdtwCrossValidation.R
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ setGeneric("twdtwCrossValidation",
#' @inheritParams twdtwCrossValidation-class
#' @aliases twdtwCrossValidation
#'
#' @describeIn
#' @describeIn twdtwCrossValidation
#' Splits the set of time series into training and validation.
#' The function uses stratified sampling and a simple random sampling for
#' each stratum. For each data partition this function performs a TWDTW
Expand Down
75 changes: 75 additions & 0 deletions R/xtable.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,75 @@
#' @aliases xtable
#' @inheritParams twdtwAssessment-class
#' @rdname twdtwAssessment-class
#'
#' @param x an object of class \code{\link[dtwSat]{twdtwAssessment}}.
#' @param type table type, 'accuracy' for User's and Producer's Accuracy,
#' 'matrix' for error matrix, and 'area' for area and uncertainty.
#' Default is 'accuracy'.
#' @param time.labels a character or numeric for the time period or NULL to
#' include all classified periods. Default is NULL.
#' @param ... other arguments to pass to \code{\link[xtable]{xtable}}.
#'
#' @export
setMethod("xtable",
signature = signature(x = "twdtwAssessment"),
definition = function(x, type="accuracy", time.labels=NULL, ...){
xtable(x@accuracySummary$ErrorMatrix, ...)
pt = pmatch(type,c("accuracy","matrix","area"))
switch(pt,
.xtable.accuracy(x=x@accuracySummary, ...),
.xtable.matrix(x=x@accuracySummary$ErrorMatrix, ...),
.xtable.area(x=x@accuracySummary$AreaUncertainty, ...)
)
}
)

.xtable.accuracy = function(x, ...){

prop = x$ProportionMatrix
prop = data.frame(apply(prop[,!names(prop)%in%c("Area","w")], 1, FUN = sprintf, fmt="%.2f"), stringsAsFactors = FALSE)
prop$`User's` = ""
prop$`Producers's` = ""
prop$`Overall` = ""

ua = sprintf("%.2f$\\pm$%.2f", round(x$UsersAccuracy[,"Accuracy"],2), round(x$UsersAccuracy[,"ci"], 2))
pa = sprintf("%.2f$\\pm$%.2f", round(x$ProducersAccuracy[,"Accuracy"],2), round(x$ProducersAccuracy[,"ci"], 2))
oa = sprintf("%.2f$\\pm$%.2f", round(x$OverallAccuracy["Accuracy"],2), round(x$OverallAccuracy["ci"], 2))

prop$`User's`[1:length(ua)] = ua
prop$`Producers's`[1:length(pa)] = pa
prop$`Overall`[1:length(oa)] = oa
tbl = xtable(prop)

comment = list()
comment$pos = list()
comment$pos[[1]] = c(nrow(tbl))
comment$command = c(paste("\\hline \n", "Test foot note. \n", sep = ""))

print.xtable(tbl, add.to.row = comment, sanitize.text.function = function(x) x)
}

.xtable.matrix = function(x, ...){
tbl = xtable(x, digits = c(rep(0, ncol(x)-1), 2, 2), ...)
print.xtable(tbl, sanitize.text.function = function(x) x)
}

.xtable.area = function(x, ...){

x=twdtw_assess@accuracySummary$AreaUncertainty
x = data.frame(x)

mp = sprintf("%.2f", round(unlist(x$Mapped),2))
ad = sprintf("%.2f", round(unlist(x$Adjusted),2))
ci = sprintf("$\\pm$%.2f", round(unlist(x$ci),2))

tbl = data.frame(mp, ad, ci)
names(tbl) = c("Mapped", "Adjusted", "Margin of error (95\\% CI)")

tbl = xtable(tbl)

print.xtable(tbl, sanitize.text.function = function(x) x)
}



1 change: 1 addition & 0 deletions R/zzz.R
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@
#' @importFrom stats xtabs ave window na.omit sd qnorm
#' @importFrom lubridate month month<- day day<- year year<-
#' @importFrom caret createDataPartition
#' @importFrom xtable xtable print.xtable
#'
NULL

Expand Down
31 changes: 31 additions & 0 deletions man/MOD13Q1.MT.yearly.patterns.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

37 changes: 37 additions & 0 deletions man/MOD13Q1.patterns.list.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

34 changes: 34 additions & 0 deletions man/MOD13Q1.ts.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

23 changes: 23 additions & 0 deletions man/MOD13Q1.ts.labels.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

Loading

0 comments on commit 59ece52

Please sign in to comment.