R/CIA.r

In intensity.analysis: Intensity of Change for Comparing Categorical Maps from Sequential Intervals

#' Category level intensity analysis
#' @details Gets the list of crosstabulation tables, time points and categories vectors and returns a list of gain and loss metrics accompanied with relevant bar graphs.
#' @param crosstabulation List of crosstabulation tables generated by \code{multicrosstab} function.
#' @param time.points a charachter vector showing the time point of each raster layer in chronological order.
#' @param categories A charachter vector showing the categories in the map. Order of categories decided bases on the equivalent IDs in the raster attribute table.
#' @return The output is a list of lists. Elements of the list include: gross.loss, gross.gain, loss.intensity, gain.intensity, Uniform.category.intensity, loss.behavior and gain.behavior.
#' @import ggplot2
#' @importFrom reshape2 melt
#' @importFrom graphics plot
#' @importFrom stats na.omit
#' @export
#' @examples
#' raster_2005 <- raster::raster(system.file("external/RASTER_2005.RST", package="intensity.analysis"))
#' raster_2010 <- raster::raster(system.file("external/RASTER_2010.RST", package="intensity.analysis"))
#' raster_2012 <- raster::raster(system.file("external/RASTER_2012.RST", package="intensity.analysis"))
#' raster.layers <- list(raster_2005, raster_2010, raster_2012)
#' time.points <- c("2005","2010","2012")
#' categories <- c("Water","Trees","Impervious")
#' crosstabulation <- multicrosstab(raster.layers, time.points, categories)
#' CIA.output <- CIA(crosstabulation, time.points, categories)



CIA <- function(crosstabulation, time.points, categories){

  # call internal function
  parameters <- reqpar(time.points)

  # Each element of the list is a placeholder of the equivalent variable for each time interval.
  gross.loss <- list()
  gross.gain <- list()
  loss.intensity <- list()
  gain.intensity <- list()
  loss.behavior <- list()
  gain.behavior <- list()

  # Uniform intensity on the category level. (This is the same as intensity in the interval level)
  Uniform.category.intensity <- vector(mode = "numeric", length = as.integer(parameters$number.of.intervals))

  Epsilon <- 0.0000001
  number.of.intervals <- as.integer(parameters$number.of.intervals)

  for (i in 1: number.of.intervals){

    # Diagonal elements of the crosstab matrix
    persistance <- diag(crosstabulation[[i]])

    # Sum of all the rows in each crosstabulation matrix (what we have at the beginning of each interval)
    sum.initial <- rowSums(crosstabulation[[i]])

    # Sum of all the columns in each crosstabulation matrix (what we have at the end of each interval)
    sum.subsequent <- colSums(crosstabulation[[i]])

    # Sum of non-diagonal elements of matrix for each row considering the duration of the interval (annual loss size)
    gross.loss[[i]] <- (sum.initial - persistance) / parameters$duration[[i]]

    # Gross loss of each category considering the size of the category at initial time (annual loss intensity).
    loss.intensity[[i]] <- gross.loss[[i]] / sum.initial

    # Sum of non-diagonal elements of matrix for each column considering the duration of the interval (annual gain size).
    gross.gain[[i]] <- (sum.subsequent - persistance) / parameters$duration[[i]]

    # Gross gain of each category considering the size of the category at the subsequent time (annual gain intensity).
    gain.intensity[[i]] <- gross.gain[[i]] / sum.subsequent

    # Total amount of change in each interval
    overal.difference <- (sum(crosstabulation[[i]]) - sum(persistance))

    # Total amount of change in regard to the geographical extent
    # of the study area and duration of each interval
    Uniform.category.intensity[i] <- overal.difference / (sum(crosstabulation[[i]]) * parameters$duration[[i]])

    # Comparison of loss intensity in each interval with the uniform intensity for each category.
    loss.behavior[[i]] <- ifelse(abs(loss.intensity[[i]] - Uniform.category.intensity[i]) < Epsilon, "Uniform", ifelse (loss.intensity[[i]] < Uniform.category.intensity[i], "Dormant", "Active"))

    # Comparison of gain intensity in each interval with the uniform intensity for each category.
    gain.behavior[[i]] <- ifelse(abs(gain.intensity[[i]] - Uniform.category.intensity[i]) < Epsilon, "Uniform",
                                 ifelse (gain.intensity[[i]] < Uniform.category.intensity[i], "Dormant", "Active"))

    ## Annual Change preparation for ggplot2 ---------------------------------
    # Make a dataframe of categories, gross loss and gross gain.
    gross <- as.data.frame(cbind(categories,gross.loss[[i]], gross.gain[[i]]))
    colnames(gross) <- c("categories","gross.loss","gross.gain")
    gross$gross.loss <- as.numeric(as.character(gross$gross.loss))
    gross$gross.gain <- as.numeric(as.character(gross$gross.gain))
    # Use the melt function from reshape2 package to plot both gross gain and loss in a single bar graph.
    gross.m <- melt(gross, id.vars="categories")

    # Force loading of "Times" font family to avoid "font family not found in Windows font database" warning message while plotting.
    # windowsFonts("Times" = windowsFont("TT Times New Roman"))


    # plot of Annual Change in number of elements unit
    plot(ggplot(gross.m, aes(x = gross.m$categories, y = gross.m$value, fill = gross.m$variable)) +
           geom_bar(width = 0.5, position = position_dodge(width=0.6), stat="identity") + theme_bw() + coord_flip() +
           labs(title= paste("Annual Change Size by Category during ", as.character(parameters$initial.times[[i]]),as.character("-"),
                             as.character(parameters$subsequent.times[[i]])), x= "Category", y= "Annual Change (# of elements)" ) +
           scale_fill_manual(labels = c("Annual Loss", "Annual Gain"), values = c("#950101", "#007a01")) +
           theme(plot.title = element_text(family = "Times", color="#353535", face="bold", size=14, hjust=0.5)) +
           theme(legend.position="bottom", legend.title = element_blank()))

    ## Annual Change intensity preparation for ggplot2-----------------------------------
    # Make a dataframe of categories, loss and gross gain intensities.
    gross.intensity <- as.data.frame(cbind(categories,loss.intensity[[i]], gain.intensity[[i]]))
    colnames(gross.intensity) <- c("categories","loss.intensity", "gain.intensity")
    gross.intensity$loss.intensity <- as.numeric(as.character(gross.intensity$loss.intensity))
    gross.intensity$gain.intensity <- as.numeric(as.character(gross.intensity$gain.intensity))
    # Use the melt function from reshape2 package to plot both gain and loss intensities in a single bar graph.
    gross.intensity.m <- melt(gross.intensity, id.vars="categories")

    # plot of Annual Change intensity in percentage of category
    plot(ggplot(gross.intensity.m, aes(x = gross.intensity.m$categories, y = gross.intensity.m$value * 100, fill = gross.intensity.m$variable)) +
           geom_bar(width = 0.5, position = position_dodge(width=0.6), stat="identity") + theme_bw() + coord_flip() +
           labs(title = paste("Annual Change Intensity by Category during ", as.character(parameters$initial.times[[i]]),as.character("-"),
                              as.character(parameters$subsequent.times[[i]])), x = "Category", y = "Annual Change Intensity (% of category)") +
           scale_fill_manual(labels = c("Annual Loss / Initial Size", "Annual Gain / Final Size"), values = c("#fe0101", "#00f801")) +
           geom_hline(aes(yintercept = Uniform.category.intensity[[i]] * 100), linetype= "dashed", colour = "blue", size = 1, show.legend = FALSE) + coord_flip() +
           geom_text(aes(x = (length(categories) / 2), y = Uniform.category.intensity[[i]] * 100, label = "Active"), family = "Times", colour = "black", angle = 90, vjust = 1.2, show.legend = NA) +
           geom_text(aes(x = (length(categories) / 2), y = Uniform.category.intensity[[i]] * 100, label = "Dormant"), family = "Times", colour = "black", angle = 90, vjust = -1, show.legend = NA) +
           theme(plot.title = element_text(family = "Times", color = "#353535", face = "bold", size = 14, hjust = 0.5)) +
           theme(legend.position = "bottom", legend.title = element_blank()))
  }

  return(list(Annual.Gross.Loss = gross.loss, Annual.Gross.Gain = gross.gain, Annual.Loss.Intensity = loss.intensity, Annual.Gain.Intensity = gain.intensity,
              Uniform.Category.Intensity = Uniform.category.intensity, Loss.Behavior = loss.behavior, Gain.Behavior = gain.behavior))
}