From a51c0349415a9ef9dd93edeb055bc08fc52a5eb5 Mon Sep 17 00:00:00 2001 From: vwmaus Date: Thu, 9 Feb 2017 16:18:28 +0100 Subject: [PATCH] Include plot methods for twdtwAssessment and twdtwCrossValidation --- DESCRIPTION | 3 +- NAMESPACE | 3 +- R/plot.R | 16 +++- R/plotAccuracy.R | 161 +++++++++++++++++++++++++++++++++++++ R/plotAdjustedArea.R | 89 ++++++++++++++++++++ R/plotCrossValidation.R | 89 -------------------- man/plot.Rd | 4 + man/plotAccuracy.Rd | 55 +++++++++++++ man/plotAdjustedArea.Rd | 51 ++++++++++++ man/plotCrossValidation.Rd | 66 --------------- 10 files changed, 379 insertions(+), 158 deletions(-) create mode 100644 R/plotAccuracy.R create mode 100644 R/plotAdjustedArea.R delete mode 100644 R/plotCrossValidation.R create mode 100644 man/plotAccuracy.Rd create mode 100644 man/plotAdjustedArea.Rd delete mode 100644 man/plotCrossValidation.Rd diff --git a/DESCRIPTION b/DESCRIPTION index ebdc21f..ecd139b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -71,12 +71,13 @@ Collate: 'methods.R' 'miscellaneous.R' 'plot.R' + 'plotAccuracy.R' + 'plotAdjustedArea.R' 'plotAlignments.R' 'plotArea.R' 'plotChanges.R' 'plotClassification.R' 'plotCostMatrix.R' - 'plotCrossValidation.R' 'plotDistance.R' 'plotMaps.R' 'plotMatches.R' diff --git a/NAMESPACE b/NAMESPACE index 755195b..1f86db4 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -6,12 +6,13 @@ export(getDatesFromDOY) export(linearWeight) export(logisticWeight) export(plot) +export(plotAccuracy) +export(plotAdjustedArea) export(plotAlignments) export(plotArea) export(plotChanges) export(plotClassification) export(plotCostMatrix) -export(plotCrossValidation) export(plotDistance) export(plotMaps) export(plotMatches) diff --git a/R/plot.R b/R/plot.R index 16d2f2c..6abd78a 100644 --- a/R/plot.R +++ b/R/plot.R @@ -54,6 +54,20 @@ #' @export NULL +#' @aliases plot-twdtwAssessment +#' @inheritParams plot +#' @rdname plot +#' @export +setMethod("plot", + signature(x = "twdtwAssessment"), + function(x, type="area", ...){ + pt = pmatch(type, c("area","accuracy")) + switch(pt, + plotAdjustedArea(x, ...), + plotAccuracy(x, ...) + ) + } +) #' @aliases plot-twdtwTimeSeries #' @inheritParams plot @@ -64,7 +78,7 @@ setMethod("plot", function(x, type="crossvalidation", ...){ pt = pmatch(type, c("crossvalidation")) switch(pt, - plotCrossValidation(x, ...) + plotAccuracy(x, ...) ) } ) diff --git a/R/plotAccuracy.R b/R/plotAccuracy.R new file mode 100644 index 0000000..ea2e7d1 --- /dev/null +++ b/R/plotAccuracy.R @@ -0,0 +1,161 @@ +#' @title Plotting accuracy assessment +#' @author Victor Maus, \email{vwmaus1@@gmail.com} +#' +#' @description Method for plotting accuracy assessment results. +#' +#' @param x An object of class \code{\link[dtwSat]{twdtwAssessment}} or +#' \code{\link[dtwSat]{twdtwCrossValidation}}. +#' +#' @param time.labels a character or numeric for the time periods or NULL to +#' aggregate all classified periods in the same plot. Default is NULL. Used +#' if \code{x} is \code{\link[dtwSat]{twdtwAssessment}}. +#' +#' @param perc if TRUE shows the results in percent of area. Otherwise shows the +#' area in the map units or km2 for no project raster. Default is TRUE. +#' +#' @param category.name a character vector defining the class names. If NULL +#' then use the classe names in the object \code{x}. Default is NULL. +#' +#' @param category.type a character defining the categories type "numeric" +#' or "letter", if NULL then use the class names. Default is NULL. +#' +#' @param conf.int confidence level (0-1) for interval estimation of the population mean. +#' for details see \code{\link[Hmisc]{smean.cl.normal}}. Used if \code{x} is +#' \code{\link[dtwSat]{twdtwCrossValidation}}. +#' +#' @return A \link[ggplot2]{ggplot} object. +#' +#' @seealso +#' \code{\link[dtwSat]{twdtwAssessment}} and \code{\link[dtwSat]{twdtwAssess}} +#' +#' @examples +#' \dontrun{ +#' +#' # See ?twdtwAssess and ?twdtwCrosValidate +#' +#' plotAccuracy(x) +#' +#' plotAccuracy(x, category.type="letter") +#' +#' } +#' +#' @export +plotAccuracy = function(x, perc=TRUE, conf.int=.95, time.labels=NULL, + category.name=NULL, category.type=NULL){ + + if(class(x)=="twdtwCrossValidation"){ + gp = .plotCrossValidation(x, conf.int, perc, category.name, category.type) + } else { + if(class(x)=="twdtwAssessment"){ + gp = .plotAssessmentAccuracy(x, perc, time.labels, category.name, category.type) + } else { + stop("class of x is not twdtwAssessment or twdtwCrossValidation") + } + } + + gp + +} + +.plotAssessmentAccuracy = function(x, perc=TRUE, time.labels=NULL, + category.name=NULL, category.type=NULL){ + + if(is.null(category.name)){ + category.name = rownames(x@accuracySummary$ProportionMatrix) + category.name = category.name[-length(category.name)] + } + if(!is.null(category.type)) + category.name = switch(pmatch(category.type,c("numeric","letter")), + as.character(seq(1:length(category.name))), + LETTERS[1:length(category.name)] + ) + + y = list(`Accumulated` = x@accuracySummary) + if(!is.null(time.labels)) + y = x@accuracyByPeriod[time.labels] + if(is.null(y)) + stop("time.labels out of bounds", call. = TRUE) + + time.labels = names(y) + if(is.null(time.labels)) + time.labels = seq_along(y) + + df = do.call("rbind", lapply(time.labels, function(i){ + # User's + df = data.frame(y[[i]]$UsersAccuracy) + df_u = data.frame(value = df$Accuracy, + variable = category.name, + Legend = "User's", + ci = df$ci, + Period = i) + # Producer's + df = data.frame(y[[i]]$ProducersAccuracy) + df_p = data.frame(value = df$Accuracy, + variable = category.name, + Legend = "Producer's", + ci = df$ci, + Period = i) + + df = rbind(df_u, df_p) + df + })) + + limits = aes_string(ymax = "value + ci", ymin = "value - ci") + dodge = position_dodge(width=0.9) + + gp = ggplot(df, aes_string(fill="Legend", y="value", x="variable")) + + facet_wrap(~Period, scales = "free") + + geom_bar(position="dodge", stat="identity", na.rm=TRUE) + + geom_errorbar(limits, position=dodge, width=0.25, na.rm=TRUE) + + scale_fill_grey(start = .6, end = .3) + + xlab("") + + ylab("Area") + + if(perc) + gp = gp + scale_y_continuous(labels = percent) + + gp + +} + +.plotCrossValidation = function(x, conf.int, perc, category.name, category.type){ + + if(is.null(category.name)){ + category.name = rownames(x@accuracy$Resample1$ErrorMatrix) + category.name = category.name + } + if(!is.null(category.type)) + category.name = switch(pmatch(category.type,c("numeric","letter")), + as.character(seq(1:length(category.name))), + LETTERS[1:length(category.name)] + ) + + UA = do.call("rbind", lapply(x@accuracy, function(x) data.frame(label="UA", rbind(x$UsersAccuracy)))) + names(UA)[-1] = category.name + PA = do.call("rbind", lapply(x@accuracy, function(x) data.frame(label="PA", rbind(x$UsersAccuracy)))) + names(PA)[-1] = category.name + df = melt(rbind(UA,PA), id="label") + df$label = factor(df$label, levels = c("UA", "PA"), + labels = c("User's Accuracy", "Producer's Accuracy")) + df$variable = factor(df$variable, levels = levels(df$variable), + labels = gsub("[.]","-",levels(df$variable))) + + gp = ggplot(df, aes_string(x="variable", y="value")) + + stat_summary(fun.data="mean_cl_boot", fun.args=list(conf.int = conf.int), + width=0.5, geom="crossbar", size=0.1, fill = "gray") + + geom_point(size=0.2) + + facet_grid(. ~ label) + + xlab("") + + ylab("Accuracy") + + coord_flip() + + if(perc){ + gp = gp + scale_y_continuous(limits = c(0,1), breaks = seq(0,1,.2), labels = percent) + } else { + gp = gp + scale_y_continuous(limits = c(0,1), breaks = seq(0,1,.2)) + } + + + gp +} + diff --git a/R/plotAdjustedArea.R b/R/plotAdjustedArea.R new file mode 100644 index 0000000..27d2ef8 --- /dev/null +++ b/R/plotAdjustedArea.R @@ -0,0 +1,89 @@ +#' @title Plotting area and uncertainty +#' @author Victor Maus, \email{vwmaus1@@gmail.com} +#' +#' @description Method for plotting area and uncertainty. +#' +#' @inheritParams plotAccuracy +#' +#' @return A \link[ggplot2]{ggplot} object. +#' +#' @seealso +#' \code{\link[dtwSat]{twdtwAssessment}} and \code{\link[dtwSat]{twdtwAssess}} +#' +#' @examples +#' \dontrun{ +#' +#' # See ?twdtwAssess +#' +#' plotAdjustedArea(twdtw_assess) +#' +#' plotAdjustedArea(twdtw_assess, category.type="letter") +#' +#' } +#' +#' @export +plotAdjustedArea = function(x, perc=TRUE, time.labels=NULL, + category.name=NULL, category.type=NULL){ + + y = list(`Accumulated area` = x@accuracySummary) + if(!is.null(time.labels)) + y = x@accuracyByPeriod[time.labels] + if(is.null(y)) + stop("time.labels out of bounds", call. = TRUE) + + if(is.null(category.name)){ + category.name = rownames(x@accuracySummary$ProportionMatrix) + category.name = category.name[-length(category.name)] + } + if(!is.null(category.type)) + category.name = switch(pmatch(category.type,c("numeric","letter")), + as.character(seq(1:length(category.name))), + LETTERS[1:length(category.name)] + ) + + time.labels = names(y) + if(is.null(time.labels)) + time.labels = seq_along(y) + + df = do.call("rbind", lapply(time.labels, function(i){ + df = data.frame(y[[i]]$AreaUncertainty) + total_area = sum(unlist(df$Mapped)) + df_m = data.frame(df$Mapped) + names(df_m) = category.name + df_m = melt(df_m) + df_m$Legend = "Mapped" + df_m$ci = as.numeric(NA) + df_m$Period = i + df_a = data.frame(df$Adjusted) + names(df_a) = category.name + df_a = melt(df_a) + df_a$Legend = "Adjusted" + df_a$ci = as.numeric(df$ci) + df_a$Period = i + df = rbind(df_m, df_a) + if(perc){ + df$value = df$value/total_area + df$ci = df$ci/total_area + } + df + })) + + limits = aes_string(ymax = "value + ci", ymin = "value - ci") + dodge = position_dodge(width=0.9) + + gp = ggplot(df, aes_string(fill="Legend", y="value", x="variable")) + + facet_wrap(~Period, scales = "free") + + geom_bar(position="dodge", stat="identity", na.rm=TRUE) + + geom_errorbar(limits, position=dodge, width=0.25, na.rm=TRUE) + + scale_fill_grey(start = .6, end = .3) + + xlab("") + + ylab("Area") + + if(perc) + gp = gp + scale_y_continuous(labels = percent) + + gp + +} + + diff --git a/R/plotCrossValidation.R b/R/plotCrossValidation.R deleted file mode 100644 index ccb1e7c..0000000 --- a/R/plotCrossValidation.R +++ /dev/null @@ -1,89 +0,0 @@ -############################################################### -# # -# (c) Victor Maus # -# Institute for Geoinformatics (IFGI) # -# University of Muenster (WWU), Germany # -# # -# Earth System Science Center (CCST) # -# National Institute for Space Research (INPE), Brazil # -# # -# # -# R Package dtwSat - 2016-11-29 # -# # -############################################################### - - -#' @title Plotting cross-validation -#' @author Victor Maus, \email{vwmaus1@@gmail.com} -#' -#' @description Method for plotting cross-validation results. -#' -#' @param x An object of class \code{\link[dtwSat]{plotCrossValidation}}. -#' @param conf.int confidence level (0-1) for interval estimation of the population mean. -#' for details see \code{\link[Hmisc]{smean.cl.normal}}. -#' -#' @return A \link[ggplot2]{ggplot} object. -#' -#' @seealso -#' \code{\link[dtwSat]{twdtwCrossValidation}} -#' -#' @examples -#' \dontrun{ -#' # Data folder -#' data_folder = system.file("lucc_MT/data", package = "dtwSat") -#' -#' # Read dates -#' dates = scan(paste(data_folder,"timeline", sep = "/"), what = "dates") -#' -#' # Read raster time series -#' evi = brick(paste(data_folder,"evi.tif", sep = "/")) -#' raster_timeseries = twdtwRaster(evi, timeline = dates) -#' -#' # Read field samples -#' field_samples = read.csv(paste(data_folder,"samples.csv", sep = "/")) -#' table(field_samples[["label"]]) -#' -#' # Read field samples projection -#' proj_str = scan(paste(data_folder,"samples_projection", sep = "/"), -#' what = "character") -#' -#' # Get sample time series from raster time series -#' field_samples_ts = getTimeSeries(raster_timeseries, -#' y = field_samples, proj4string = proj_str) -#' field_samples_ts -#' -#' # Run cross validation -#' set.seed(1) -#' # Define TWDTW weight function -#' log_fun = logisticWeight(alpha=-0.1, beta=50) -#' cross_validation = twdtwCrossValidation(field_samples_ts, times=3, p=0.1, -#' freq = 8, formula = y ~ s(x, bs="cc"), weight.fun = log_fun) -#' -#' summary(cross_validation) -#' -#' plot(cross_validation, conf.int=.99) -#' -#' } -#' -#' @export -plotCrossValidation = function(x, conf.int=.95){ - UA = do.call("rbind", lapply(x@accuracy, function(x) data.frame(label="UA", rbind(x$UsersAccuracy)))) - PA = do.call("rbind", lapply(x@accuracy, function(x) data.frame(label="PA", rbind(x$UsersAccuracy)))) - df = melt(rbind(UA,PA), id="label") - df$label = factor(df$label, levels = c("UA", "PA"), - labels = c("User's Accuracy", "Producer's Accuracy")) - df$variable = factor(df$variable, levels = levels(df$variable), - labels = gsub("[.]","-",levels(df$variable))) - gp = ggplot(df, aes_string(x="variable", y="value")) + - stat_summary(fun.data="mean_cl_boot", fun.args=list(conf.int = conf.int), - width=0.5, geom="crossbar", size=0.1, fill = "gray") + - geom_point(size=0.2) + - facet_grid(. ~ label) + - scale_y_continuous(limits = c(0,1), labels = percent, breaks = seq(0,1,.2)) + - xlab("") + - ylab("Accuracy") + - coord_flip() - gp -} - - diff --git a/man/plot.Rd b/man/plot.Rd index 49c347f..5261935 100644 --- a/man/plot.Rd +++ b/man/plot.Rd @@ -3,15 +3,19 @@ \docType{methods} \name{plot} \alias{plot} +\alias{plot,twdtwAssessment,ANY-method} \alias{plot,twdtwCrossValidation,ANY-method} \alias{plot,twdtwMatches,ANY-method} \alias{plot,twdtwRaster,ANY-method} \alias{plot,twdtwTimeSeries,ANY-method} +\alias{plot-twdtwAssessment} \alias{plot-twdtwMatches} \alias{plot-twdtwRaster} \alias{plot-twdtwTimeSeries} \title{Plotting twdtw* objects} \usage{ +\S4method{plot}{twdtwAssessment,ANY}(x, type = "area", ...) + \S4method{plot}{twdtwCrossValidation,ANY}(x, type = "crossvalidation", ...) \S4method{plot}{twdtwTimeSeries,ANY}(x, type = "timeseries", ...) diff --git a/man/plotAccuracy.Rd b/man/plotAccuracy.Rd new file mode 100644 index 0000000..1bf8dca --- /dev/null +++ b/man/plotAccuracy.Rd @@ -0,0 +1,55 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plotAccuracy.R +\name{plotAccuracy} +\alias{plotAccuracy} +\title{Plotting accuracy assessment} +\usage{ +plotAccuracy(x, perc = TRUE, conf.int = 0.95, time.labels = NULL, + category.name = NULL, category.type = NULL) +} +\arguments{ +\item{x}{An object of class \code{\link[dtwSat]{twdtwAssessment}} or +\code{\link[dtwSat]{twdtwCrossValidation}}.} + +\item{perc}{if TRUE shows the results in percent of area. Otherwise shows the +area in the map units or km2 for no project raster. Default is TRUE.} + +\item{conf.int}{confidence level (0-1) for interval estimation of the population mean. +for details see \code{\link[Hmisc]{smean.cl.normal}}. Used if \code{x} is +\code{\link[dtwSat]{twdtwCrossValidation}}.} + +\item{time.labels}{a character or numeric for the time periods or NULL to +aggregate all classified periods in the same plot. Default is NULL. Used +if \code{x} is \code{\link[dtwSat]{twdtwAssessment}}.} + +\item{category.name}{a character vector defining the class names. If NULL +then use the classe names in the object \code{x}. Default is NULL.} + +\item{category.type}{a character defining the categories type "numeric" +or "letter", if NULL then use the class names. Default is NULL.} +} +\value{ +A \link[ggplot2]{ggplot} object. +} +\description{ +Method for plotting accuracy assessment results. +} +\examples{ +\dontrun{ + +# See ?twdtwAssess and ?twdtwCrosValidate + +plotAccuracy(x) + +plotAccuracy(x, category.type="letter") + +} + +} +\author{ +Victor Maus, \email{vwmaus1@gmail.com} +} +\seealso{ +\code{\link[dtwSat]{twdtwAssessment}} and \code{\link[dtwSat]{twdtwAssess}} +} + diff --git a/man/plotAdjustedArea.Rd b/man/plotAdjustedArea.Rd new file mode 100644 index 0000000..0e0a534 --- /dev/null +++ b/man/plotAdjustedArea.Rd @@ -0,0 +1,51 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plotAdjustedArea.R +\name{plotAdjustedArea} +\alias{plotAdjustedArea} +\title{Plotting area and uncertainty} +\usage{ +plotAdjustedArea(x, perc = TRUE, time.labels = NULL, category.name = NULL, + category.type = NULL) +} +\arguments{ +\item{x}{An object of class \code{\link[dtwSat]{twdtwAssessment}} or +\code{\link[dtwSat]{twdtwCrossValidation}}.} + +\item{perc}{if TRUE shows the results in percent of area. Otherwise shows the +area in the map units or km2 for no project raster. Default is TRUE.} + +\item{time.labels}{a character or numeric for the time periods or NULL to +aggregate all classified periods in the same plot. Default is NULL. Used +if \code{x} is \code{\link[dtwSat]{twdtwAssessment}}.} + +\item{category.name}{a character vector defining the class names. If NULL +then use the classe names in the object \code{x}. Default is NULL.} + +\item{category.type}{a character defining the categories type "numeric" +or "letter", if NULL then use the class names. Default is NULL.} +} +\value{ +A \link[ggplot2]{ggplot} object. +} +\description{ +Method for plotting area and uncertainty. +} +\examples{ +\dontrun{ + +# See ?twdtwAssess + +plotAdjustedArea(twdtw_assess) + +plotAdjustedArea(twdtw_assess, category.type="letter") + +} + +} +\author{ +Victor Maus, \email{vwmaus1@gmail.com} +} +\seealso{ +\code{\link[dtwSat]{twdtwAssessment}} and \code{\link[dtwSat]{twdtwAssess}} +} + diff --git a/man/plotCrossValidation.Rd b/man/plotCrossValidation.Rd deleted file mode 100644 index 67815a0..0000000 --- a/man/plotCrossValidation.Rd +++ /dev/null @@ -1,66 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/plotCrossValidation.R -\name{plotCrossValidation} -\alias{plotCrossValidation} -\title{Plotting cross-validation} -\usage{ -plotCrossValidation(x, conf.int = 0.95) -} -\arguments{ -\item{x}{An object of class \code{\link[dtwSat]{plotCrossValidation}}.} - -\item{conf.int}{confidence level (0-1) for interval estimation of the population mean. -for details see \code{\link[Hmisc]{smean.cl.normal}}.} -} -\value{ -A \link[ggplot2]{ggplot} object. -} -\description{ -Method for plotting cross-validation results. -} -\examples{ -\dontrun{ -# Data folder -data_folder = system.file("lucc_MT/data", package = "dtwSat") - -# Read dates -dates = scan(paste(data_folder,"timeline", sep = "/"), what = "dates") - -# Read raster time series -evi = brick(paste(data_folder,"evi.tif", sep = "/")) -raster_timeseries = twdtwRaster(evi, timeline = dates) - -# Read field samples -field_samples = read.csv(paste(data_folder,"samples.csv", sep = "/")) -table(field_samples[["label"]]) - -# Read field samples projection -proj_str = scan(paste(data_folder,"samples_projection", sep = "/"), - what = "character") - -# Get sample time series from raster time series -field_samples_ts = getTimeSeries(raster_timeseries, - y = field_samples, proj4string = proj_str) -field_samples_ts - -# Run cross validation -set.seed(1) -# Define TWDTW weight function -log_fun = logisticWeight(alpha=-0.1, beta=50) -cross_validation = twdtwCrossValidation(field_samples_ts, times=3, p=0.1, - freq = 8, formula = y ~ s(x, bs="cc"), weight.fun = log_fun) - -summary(cross_validation) - -plot(cross_validation, conf.int=.99) - -} - -} -\author{ -Victor Maus, \email{vwmaus1@gmail.com} -} -\seealso{ -\code{\link[dtwSat]{twdtwCrossValidation}} -} -