diff --git a/NAMESPACE b/NAMESPACE index 3272e0d..f2451ca 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -60,7 +60,7 @@ exportMethods(show) exportMethods(subset) exportMethods(summary) exportMethods(twdtwApply) -exportMethods(twdtwAssessment) +exportMethods(twdtwAssess) exportMethods(twdtwClassify) exportMethods(twdtwCrossValidation) exportMethods(twdtwMatches) @@ -105,6 +105,7 @@ importFrom(sp,over) importFrom(sp,spTransform) importFrom(stats,ave) importFrom(stats,na.omit) +importFrom(stats,qnorm) importFrom(stats,sd) importFrom(stats,window) importFrom(stats,xtabs) diff --git a/R/class-crossValidation.R b/R/class-crossValidation.R index ba9aaf2..cd28411 100644 --- a/R/class-crossValidation.R +++ b/R/class-crossValidation.R @@ -55,10 +55,10 @@ setClass( slots = c(partitions = "list", accuracy = "list"), validity = function(object){ if(!is(object@partitions, "list")){ - stop("[twdtwTimeSeries: validation] Invalid partitions, class different from list.") + stop("[twdtwCrossValidation: validation] Invalid partitions, class different from list.") }else{} if(!is(object@accuracy, "list")){ - stop("[twdtwTimeSeries: validation] Invalid accuracy, class different from list.") + stop("[twdtwCrossValidation: validation] Invalid accuracy, class different from list.") }else{} return(TRUE) } diff --git a/R/class-twdtwAccuracy.R b/R/class-twdtwAccuracy.R index e779122..f5497a4 100644 --- a/R/class-twdtwAccuracy.R +++ b/R/class-twdtwAccuracy.R @@ -33,10 +33,15 @@ #' @param id.labels a numeric or character with an column name from \code{y} to #' be used as samples labels. Optional. #' +#' @param labels character vector with time series labels. For signature +#' \code{\link[dtwSat]{twdtwRaster}} this argument can be used to set the +#' labels for each sample in \code{y}, or it can be combined with \code{id.labels} +#' to select samples with a specific label. +#' #' @param proj4string projection string, see \code{\link[sp]{CRS-class}}. Used #' if \code{y} is a \code{\link[base]{data.frame}}. #' -#' @param conf.int specifies the confidence level. +#' @param conf.int specifies the confidence level (0-1). #' #' @seealso #' \code{\link[dtwSat]{twdtwRaster-class}}, and @@ -44,9 +49,12 @@ #' #' @section Slots : #' \describe{ -#' \item{\code{accuracy}:}{A list with the accuracy for each classified time period.} -#' \item{\code{data}:}{A \code{\link[base]{data.frame}} with reference labels, predicted labels, -#' and other TWDTW information.} +#' \item{\code{accuracySummary}:}{Overall Accuracy, User's Accuracy, Produce's Accuracy, +#' and Error Matrix (confusion matrix) considering all time periods.} +#' \item{\code{accuracyByPeriod}:}{Overall Accuracy, User's Accuracy, Produce's Accuracy, +#' and Error Matrix (confusion matrix) for each time periods independently from each other.} +#' \item{\code{data}:}{A \code{\link[base]{data.frame}} with period (from - to), reference labels, +#' predicted labels, and other TWDTW information.} #' } #' #' @examples @@ -56,13 +64,16 @@ NULL setClass( Class = "twdtwAssessment", - slots = c(accuracy = "list", data = "list"), + slots = c(accuracySummary = "list", accuracyByPeriod = "list", data = "data.frame"), validity = function(object){ - if(!is(object@partitions, "list")){ - stop("[twdtwTimeSeries: validation] Invalid partitions, class different from list.") + if(!is(object@accuracySummary, "list")){ + stop("[twdtwAssessment: validation] Invalid partitions, class different from list.") + }else{} + if(!is(object@accuracyByPeriod, "list")){ + stop("[twdtwAssessment: validation] Invalid accuracy, class different from list.") }else{} - if(!is(object@accuracy, "list")){ - stop("[twdtwTimeSeries: validation] Invalid accuracy, class different from list.") + if(!is(object@data, "data.frame")){ + stop("[twdtwAssessment: validation] Invalid accuracy, class different from data.frame.") }else{} return(TRUE) } @@ -71,14 +82,17 @@ setClass( setMethod("initialize", signature = "twdtwAssessment", definition = - function(.Object, partitions, accuracy){ - .Object@partitions = list(Resample1=NULL) - .Object@accuracy = list(OverallAccuracy=NULL, UsersAccuracy=NULL, ProducersAccuracy=NULL, - error.matrix=table(NULL), data=data.frame(NULL)) - if(!missing(partitions)) - .Object@partitions = partitions - if(!missing(accuracy)) - .Object@accuracy = accuracy + function(.Object, accuracySummary, accuracyByPeriod, data){ + .Object@accuracySummary = list(OverallAccuracy=NULL, UsersAccuracy=NULL, ProducersAccuracy=NULL, ErrorMatrix=table(NULL)) + .Object@accuracyByPeriod = list(list(OverallAccuracy=NULL, UsersAccuracy=NULL, ProducersAccuracy=NULL, + ErrorMatrix=table(NULL))) + .Object@data = data.frame(Period=NULL, from=NULL, to=NULL, Distance=NULL, Predicted=NULL, Reference=NULL) + if(!missing(accuracySummary)) + .Object@accuracySummary = accuracySummary + if(!missing(accuracyByPeriod)) + .Object@accuracyByPeriod = accuracyByPeriod + if(!missing(data)) + .Object@data = data validObject(.Object) return(.Object) } diff --git a/R/getTimeSeries.R b/R/getTimeSeries.R index ef03a4b..131ac35 100644 --- a/R/getTimeSeries.R +++ b/R/getTimeSeries.R @@ -116,29 +116,15 @@ getTimeSeries.twdtwTimeSeries = function(object, labels){ setMethod("getTimeSeries", "twdtwRaster", function(object, y, labels=NULL, proj4string = NULL, id.labels=NULL){ - if(!"label"%in%names(y)) y$label = paste0("ts",row.names(y)) - if(!is.null(id.labels)) y$label = as.character(y[[id.labels]]) - if(!is.null(id.labels) & !is.null(labels)){ - I = which(!is.na(match(as.character(y$label), as.character(labels)))) - if(length(I)<1) - stop("there is no matches between id.labels and labels") - } else if(!is.null(labels)) { - y$label = as.character(labels) - } + y = .adjustLabelID(y, labels, id.labels) + if(!"from"%in%names(y)) y$from = as.Date(index(object)[1]) if(!"to"%in%names(y)) y$to = as.Date(tail(index(object),1)) - if(is(y, "data.frame")){ - if(is.null(proj4string)){ - warning("Missing projection. Setting the same projection as the raster time series.", call. = FALSE) - proj4string = CRS(projection(object)) - } - if(!is(proj4string, "CRS")) proj4string = try(CRS(proj4string)) - y = SpatialPointsDataFrame(y[,c("longitude","latitude")], y, proj4string = proj4string) - } - if(!(is(y, "SpatialPoints") | is(y, "SpatialPointsDataFrame"))) - stop("y is not SpatialPoints or SpatialPointsDataFrame") + + y = .toSpatialPointsDataFrame(y, object, proj4string) + extractTimeSeries.twdtwRaster(object, y) }) diff --git a/R/miscellaneous.R b/R/miscellaneous.R index 55b5dbe..47c9d22 100644 --- a/R/miscellaneous.R +++ b/R/miscellaneous.R @@ -117,4 +117,30 @@ shiftDates.twdtwTimeSeries = function(x, year){ data = data.frame(Predicted=pred, Reference=ref) } +.adjustLabelID = function(y, labels, id.labels){ + if(!"label"%in%names(y)) y$label = paste0("ts",row.names(y)) + if(!is.null(id.labels)) y$label = as.character(y[[id.labels]]) + if(!is.null(id.labels) & !is.null(labels)){ + I = which(!is.na(match(as.character(y$label), as.character(labels)))) + if(length(I)<1) + stop("there is no matches between id.labels and labels") + } else if(!is.null(labels)) { + y$label = as.character(labels) + } + y +} +.toSpatialPointsDataFrame = function(y, object, proj4string){ + if(is(y, "data.frame")){ + if(is.null(proj4string)){ + warning("Missing projection. Setting the same projection as the raster time series.", call. = FALSE) + proj4string = CRS(projection(object)) + } + if(!is(proj4string, "CRS")) proj4string = try(CRS(proj4string)) + y = SpatialPointsDataFrame(y[,c("longitude","latitude")], y, proj4string = proj4string) + } + if(!(is(y, "SpatialPoints") | is(y, "SpatialPointsDataFrame"))) + stop("y is not SpatialPoints or SpatialPointsDataFrame") + row.names(y) = 1:nrow(y) + y +} diff --git a/R/twdtwAssessment.R b/R/twdtwAssessment.R index 26a6023..08529b3 100644 --- a/R/twdtwAssessment.R +++ b/R/twdtwAssessment.R @@ -1,10 +1,10 @@ -setGeneric("twdtwAssessment", - def = function(object, ...) standardGeneric("twdtwAssessment") +setGeneric("twdtwAssess", + def = function(object, ...) standardGeneric("twdtwAssess") ) #' @inheritParams twdtwAssessment-class -#' @aliases twdtwAssessment +#' @aliases twdtwAssess #' #' @describeIn twdtwAssessment this function performs an accuracy assessment #' of the classified maps. The function returns Overall Accuracy, @@ -24,46 +24,59 @@ setGeneric("twdtwAssessment", #' mir = brick(system.file("lucc_MT/data/mir.tif", package="dtwSat")) #' doy = brick(system.file("lucc_MT/data/doy.tif", package="dtwSat")) #' timeline = scan(system.file("lucc_MT/data/timeline", package="dtwSat"), what="date") -#' #' rts = twdtwRaster(evi, ndvi, red, blue, nir, mir, timeline = timeline, doy = doy) +#' +#' # Read fiels samples #' field_samples = read.csv(system.file("lucc_MT/data/samples.csv", package="dtwSat")) #' proj_str = scan(system.file("lucc_MT/data/samples_projection", #' package="dtwSat"), what = "character") -#' field_samples_ts = getTimeSeries(rts, y = field_samples, proj4string = proj_str) -#' temporal_patterns = createPatterns(field_samples_ts, freq = 8, formula = y ~ s(x)) -#' log_fun = weight.fun=logisticWeight(-0.1,50) #' -#' # Run TWDTW analysis for raster time series +#' # Split samples for training (10%) and validation (90%) using stratified sampling +#' library(caret) +#' set.seed(1) +#' I = unlist(createDataPartition(field_samples$label, p = 0.1)) +#' training_samples = field_samples[I,] +#' validation_samples = field_samples[-I,] #' -#' r_twdtw = twdtwApply(x=rts, y=temporal_patterns, weight.fun=log_fun, format="GTiff", -#' overwrite=TRUE, chunk.size=1000) -#' -#' r_lucc = twdtwClassify(r_twdtw, format="GTiff", overwrite=TRUE) +#' # Create temporal patterns +#' training_ts = getTimeSeries(rts, y = training_samples, proj4string = proj_str) +#' temporal_patterns = createPatterns(training_ts, freq = 8, formula = y ~ s(x)) #' +#' # Run TWDTW analysis for raster time series #' log_fun = weight.fun=logisticWeight(-0.1,50) -#' time_interval = seq(from=as.Date("2007-09-01"), to=as.Date("2013-09-01"), -#' by="12 month") -#' r_twdtw = twdtwApply(x=rts, y=patt, weight.fun=log_fun, breaks=time_interval, -#' filepath="~/test_twdtw", overwrite=TRUE, format="GTiff") -#' -#' r_lucc = twdtwClassify(r_twdtw, format="GTiff", overwrite=TRUE) -#' -#' plotMaps(r_lucc) -#' -#' # Map assessment -#' +#' r_twdtw = twdtwApply(x=rts, y=temporal_patterns, weight.fun=log_fun, format="GTiff", +#' overwrite=TRUE) +#' +#' # Classify raster based on the TWDTW analysis +#' r_lucc = twdtwClassify(r_twdtw, format="GTiff", overwrite=TRUE, filepath="res1") +#' plot(r_lucc) #' -#' +#' # Assess classification +#' twdtw_assess = twdtwAssess(r_lucc, validation_samples, proj4string=proj_str) +#' twdtw_assess@data +#' #' } #' @export -setMethod(f = "twdtwAssessment", - definition = function(object, y, id.labels, proj4string, conf.int) - twdtwAssessment.twdtwRaster(object, y, id.labels, proj4string, conf.int)) - -twdtwAssessment.twdtwRaster = function(object, y, id.labels, proj4string, conf.int){ +setMethod(f = "twdtwAssess", signature = "twdtwRaster", + definition = function(object, y, labels=NULL, id.labels=NULL, proj4string=NULL, conf.int=.95) + twdtwAssess.twdtwRaster(object, y, labels, id.labels, proj4string, conf.int)) +twdtwAssess.twdtwRaster = function(object, y, labels, id.labels, proj4string, conf.int){ + + # Check control points + y = .adjustLabelID(y, labels, id.labels) + if(!"from"%in%names(y)) + stop("samples starting date not found, the argument 'y' must have a column called 'from'") + if(!"to"%in%names(y)) + stop("samples ending date not found, the argument 'y' must have a column called 'to'") + y = .toSpatialPointsDataFrame(y, object, proj4string) + # Get classified raster x = object@timeseries$Class + x_twdtw = object@timeseries$Distance + + # Reproject points to raster projection + y = spTransform(y, CRS(projection(x))) # Get time intervals timeline = index(object) @@ -75,7 +88,7 @@ twdtwAssessment.twdtwRaster = function(object, y, id.labels, proj4string, conf.i rlevels = levels(object) # Compute area of each class by classification interval - a_by_interval = lapply(1:nlayers(x), FUN = .area_by_class, x, rlevels, rnames) + a_by_interval = lapply(1:nlayers(x), FUN = .getAreaByClass, x, rlevels, rnames) # Compute total area by class area_by_class = do.call("rbind", a_by_interval) @@ -83,45 +96,25 @@ twdtwAssessment.twdtwRaster = function(object, y, id.labels, proj4string, conf.i # Get classified and predicted land cover/use classes for each control point pred_classes = extract(x, y) - samples_by_year = lapply(1:nrow(r_intervals), FUN = .get_pre_ref_classes, r_intervals, pred_classes, y, rlevels, rnames) - samples_all = do.call("rbind", samples_by_year) + pred_distance = extract(x_twdtw, y) + samples_by_period = lapply(1:nrow(r_intervals), FUN = .getPredRefClasses, r_intervals, pred_classes, pred_distance, y, rlevels, rnames) + samples_all = do.call("rbind", samples_by_period) # Compute error matrix - error_matrix_by_year = lapply(samples_by_year, table) - error_matrix_summary = table(samples_all) + error_matrix_by_period = lapply(1:nrow(r_intervals), function(i) table(samples_by_period[[i]][,c("Predicted","Reference")])) + error_matrix_summary = table(samples_all[,c("Predicted","Reference")]) # Compute accuracy assessment - accuracy_by_year = lapply(seq_along(error_matrix_by_year), function(i) .twdtwAssessment(x = error_matrix_by_year[[i]], a_by_interval[[i]], conf.int)) - accuracy_summary = .twdtwAssessment(error_matrix_summary, area_by_class, conf.int=1.96) + accuracy_by_period = lapply(seq_along(error_matrix_by_period), function(i) .twdtwAssess(x = error_matrix_by_period[[i]], a_by_interval[[i]], conf.int=conf.int)) + accuracy_summary = .twdtwAssess(error_matrix_summary, area_by_class, conf.int=conf.int) - # new("twdtwCrossValidation", partitions=partitions, accuracy=res) + new("twdtwAssessment", accuracySummary=accuracy_summary, accuracyByPeriod=accuracy_by_period, data=samples_all) - list(accuracy_summary, accuracy_by_year) - -} - -.get_pre_ref_classes = function(i, r_intervals, pred, y, rlevels, rnames){ - I = which((r_intervals$to[i] - as.Date(y$from) > 30) & (as.Date(y$to) - r_intervals$from[i] > 30) ) - if(length(I)<1) - return(NULL) - J = match(pred[I,i], rlevels) - Predicted = factor(as.character(rnames[J]), levels = rnames, labels = rnames) - Reference = factor(as.character(y$label[I]), levels = rnames, labels = rnames) - data.frame(Predicted, Reference) } -.area_by_class = function(l, r, rlevels, rnames){ - r = raster(r, layer = l) - a = zonal(r, r, 'count') - I = match(a[,'zone'], rlevels) - out = rep(0, length(rnames)) - names(out) = rnames - out[I] = a[,'count'] * prod(res(r)) - names(out) = rnames - out -} - -.twdtwAssessment = function(x, area, conf.int){ +.twdtwAssess = function(x, area, conf.int){ + + mult = qnorm(1-(1-conf.int)/2, mean = 0, sd = 1) cnames = names(area) # cnames = paste0("aux_classname_",seq_along(cnames)) @@ -167,16 +160,15 @@ twdtwAssessment.twdtwRaster = function(object, y, id.labels, proj4string, conf.i #a_pixel = as.numeric(prop_matrix["Total",cnames] * prop_matrix["Total","A"]) #names(a_pixel) = cnames #a_ha = a_pixel*res^2/100^2 - conf.int = 1.96 temp = w^2*UA*(1-UA)/(total_map-1) VO = sum(temp, na.rm = TRUE) SO = sqrt(VO) - OCI = SO * conf.int + OCI = SO * mult VU = UA*(1-UA)/(total_map-1) SU = sqrt(VU) - UCI = SU * conf.int + UCI = SU * mult fun1 = function(x, xt, A){ sum(A*x/xt, na.rm = TRUE) @@ -197,7 +189,7 @@ twdtwAssessment.twdtwRaster = function(object, y, id.labels, proj4string, conf.i # VP = (1/Nj^2)*(expr1+expr2) 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 * 1.96 + PCI = SP * mult res = list(OverallAccuracy = c(Accuracy=OA, Var=VO, sd=SO, ci95=OCI), UsersAccuracy = cbind(Accuracy=UA, Var=VU, sd=SU, ci95=UCI), @@ -209,6 +201,29 @@ twdtwAssessment.twdtwRaster = function(object, y, id.labels, proj4string, conf.i } +.getPredRefClasses = function(i, r_intervals, pred, pred_distance, y, rlevels, rnames){ + I = which((r_intervals$to[i] - as.Date(y$from) > 30) & (as.Date(y$to) - r_intervals$from[i] > 30) ) + if(length(I)<1) + return(NULL) + J = match(pred[I,i], rlevels) + Predicted = factor(as.character(rnames[J]), levels = rnames, labels = rnames) + Reference = factor(as.character(y$label[I]), levels = rnames, labels = rnames) + #d = pred_distance[J] + data.frame(Period=i, from=r_intervals$from[i], to=r_intervals$to[i], Predicted, Reference) +} + +.getAreaByClass = function(l, r, rlevels, rnames){ + r = raster(r, layer = l) + a = zonal(r, r, 'count') + I = match(a[,'zone'], rlevels) + out = rep(0, length(rnames)) + names(out) = rnames + out[I] = a[,'count'] * prod(res(r)) + names(out) = rnames + out +} + + diff --git a/R/zzz.R b/R/zzz.R index d571226..96a76dc 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -32,7 +32,7 @@ #' @importFrom sp Polygon Polygons SpatialPoints SpatialPolygons SpatialPointsDataFrame over CRS spTransform coordinates bbox #' @importFrom mgcv gam predict.gam #' @importFrom RColorBrewer brewer.pal -#' @importFrom stats xtabs ave window na.omit sd +#' @importFrom stats xtabs ave window na.omit sd qnorm #' @importFrom lubridate month month<- day day<- year year<- #' @importFrom caret createDataPartition #' diff --git a/dtwSat.Rproj b/dtwSat.Rproj index e4046f4..a72b742 100644 --- a/dtwSat.Rproj +++ b/dtwSat.Rproj @@ -17,5 +17,5 @@ PackageUseDevtools: Yes PackageInstallArgs: --no-multiarch --with-keep.source PackageBuildArgs: --compact-vignettes=gs+qpdf PackageBuildBinaryArgs: --compact-vignettes=gs+qpdf -PackageCheckArgs: --as-cran +PackageCheckArgs: --as-cran --no-build-vignettes --no-vignettes PackageRoxygenize: rd,collate diff --git a/man/twdtwAssessment-class.Rd b/man/twdtwAssessment-class.Rd index 7ec06a6..ca28222 100644 --- a/man/twdtwAssessment-class.Rd +++ b/man/twdtwAssessment-class.Rd @@ -2,12 +2,14 @@ % Please edit documentation in R/class-twdtwAccuracy.R, R/twdtwAssessment.R \docType{methods} \name{twdtwAssessment-class} +\alias{twdtwAssess} +\alias{twdtwAssess,twdtwRaster-method} \alias{twdtwAssessment} -\alias{twdtwAssessment,ANY-method} \alias{twdtwAssessment-class} \title{class "twdtwAssessment"} \usage{ -\S4method{twdtwAssessment}{ANY}(object, y, id.labels, proj4string, conf.int) +\S4method{twdtwAssess}{twdtwRaster}(object, y, labels = NULL, + id.labels = NULL, proj4string = NULL, conf.int = 0.95) } \arguments{ \item{object}{an object of class \code{\link[dtwSat]{twdtwRaster}} resulting from @@ -20,20 +22,25 @@ whose attributes are the start ''from'' and the end ''to'' of the time interval. If missing ''from'' and/or ''to'', they are set to the time range of the \code{object}.} +\item{labels}{character vector with time series labels. For signature +\code{\link[dtwSat]{twdtwRaster}} this argument can be used to set the +labels for each sample in \code{y}, or it can be combined with \code{id.labels} +to select samples with a specific label.} + \item{id.labels}{a numeric or character with an column name from \code{y} to be used as samples labels. Optional.} \item{proj4string}{projection string, see \code{\link[sp]{CRS-class}}. Used if \code{y} is a \code{\link[base]{data.frame}}.} -\item{conf.int}{specifies the confidence level.} +\item{conf.int}{specifies the confidence level (0-1).} } \description{ This class stores the map assessment. } \section{Methods (by generic)}{ \itemize{ -\item \code{twdtwAssessment}: this function performs an accuracy assessment +\item \code{twdtwAssess}: this function performs an accuracy assessment of the classified maps. The function returns Overall Accuracy, User's Accuracy, Produce's Accuracy, and error matrix (confusion matrix) for each time interval and a summary considering all classified intervals. @@ -41,9 +48,12 @@ each time interval and a summary considering all classified intervals. \section{Slots }{ \describe{ - \item{\code{accuracy}:}{A list with the accuracy for each classified time period.} - \item{\code{data}:}{A \code{\link[base]{data.frame}} with reference labels, predicted labels, - and other TWDTW information.} + \item{\code{accuracySummary}:}{Overall Accuracy, User's Accuracy, Produce's Accuracy, + and Error Matrix (confusion matrix) considering all time periods.} + \item{\code{accuracyByPeriod}:}{Overall Accuracy, User's Accuracy, Produce's Accuracy, + and Error Matrix (confusion matrix) for each time periods independently from each other.} + \item{\code{data}:}{A \code{\link[base]{data.frame}} with period (from - to), reference labels, + predicted labels, and other TWDTW information.} } } \examples{ @@ -62,36 +72,37 @@ nir = brick(system.file("lucc_MT/data/nir.tif", package="dtwSat")) mir = brick(system.file("lucc_MT/data/mir.tif", package="dtwSat")) doy = brick(system.file("lucc_MT/data/doy.tif", package="dtwSat")) timeline = scan(system.file("lucc_MT/data/timeline", package="dtwSat"), what="date") - rts = twdtwRaster(evi, ndvi, red, blue, nir, mir, timeline = timeline, doy = doy) + +# Read fiels samples field_samples = read.csv(system.file("lucc_MT/data/samples.csv", package="dtwSat")) proj_str = scan(system.file("lucc_MT/data/samples_projection", package="dtwSat"), what = "character") -field_samples_ts = getTimeSeries(rts, y = field_samples, proj4string = proj_str) -temporal_patterns = createPatterns(field_samples_ts, freq = 8, formula = y ~ s(x)) -log_fun = weight.fun=logisticWeight(-0.1,50) -# Run TWDTW analysis for raster time series - -r_twdtw = twdtwApply(x=rts, y=temporal_patterns, weight.fun=log_fun, format="GTiff", - overwrite=TRUE, chunk.size=1000) +# Split samples for training (10\%) and validation (90\%) using stratified sampling +library(caret) +set.seed(1) +I = unlist(createDataPartition(field_samples$label, p = 0.1)) +training_samples = field_samples[I,] +validation_samples = field_samples[-I,] -r_lucc = twdtwClassify(r_twdtw, format="GTiff", overwrite=TRUE) +# Create temporal patterns +training_ts = getTimeSeries(rts, y = training_samples, proj4string = proj_str) +temporal_patterns = createPatterns(training_ts, freq = 8, formula = y ~ s(x)) +# Run TWDTW analysis for raster time series log_fun = weight.fun=logisticWeight(-0.1,50) -time_interval = seq(from=as.Date("2007-09-01"), to=as.Date("2013-09-01"), - by="12 month") -r_twdtw = twdtwApply(x=rts, y=patt, weight.fun=log_fun, breaks=time_interval, - filepath="~/test_twdtw", overwrite=TRUE, format="GTiff") - -r_lucc = twdtwClassify(r_twdtw, format="GTiff", overwrite=TRUE) - -plotMaps(r_lucc) - -# Map assessment - - - +r_twdtw = twdtwApply(x=rts, y=temporal_patterns, weight.fun=log_fun, format="GTiff", + overwrite=TRUE) + +# Classify raster based on the TWDTW analysis +r_lucc = twdtwClassify(r_twdtw, format="GTiff", overwrite=TRUE, filepath="res1") +plot(r_lucc) + +# Assess classification +twdtw_assess = twdtwAssess(r_lucc, validation_samples, proj4string=proj_str) +twdtw_assess@data + } } \author{