R/TIA.r

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

#' Transition 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: transition intensity, uniform transition, transition behavior for gain of a category. These metrics are calculated for each interval.
#' @import ggplot2
#' @importFrom stats na.omit
#' @importFrom graphics plot
#' @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)
#' TIA.output <- TIA(crosstabulation, time.points, categories)



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

  parameters <- reqpar(time.points)

  # Each element of the list is a placeholder of the equivalent variable for each time interval.
  transition.intensity <- list()
  uniform.transition <- list()
  Annual.transition.size <- list()
  transition.behavior.for.gain <- list()
  Epsilon <- 0.0000001
  number.of.intervals <- as.integer(parameters$number.of.intervals)


  for (i in 1: number.of.intervals){
    # Calculate the transitin intensity matrix by dividing the crosstab matrix by sum of the
    # categories at initial time (sum of rows) with respect to the duration of interval.
    transition.intensity[[i]] <- crosstabulation[[i]] / (rowSums(crosstabulation[[i]]) * parameters$duration[i])

    # replace the diagonal of the matrix with "NA".
    diag(transition.intensity[[i]]) <- NA

    # Uniform transitin intensity is calculated by dividing the Gain of a category to
    # the number of pixels that are not belong to that category at initial time.
    uniform.transition[[i]] <- ((colSums(crosstabulation[[i]]) - diag(crosstabulation[[i]])) / parameters$duration[[i]]) / (sum(crosstabulation[[i]]) - rowSums(crosstabulation[[i]]))

    # Comparing the transitin intensity for gain of a category and uniform transitin intensity
    # of that category to find if the Gain for that specific category is uniform, targeting or
    # avoiding the other categories loss.
    # If transitin intensity for gain of a category is equal to the uniform trnsitin intensity it labels as "uniform".
    # If transitin intensity for gain of a category is bigger in value the uniform transitin intensity it labels as "Targeting".
    # If transitin intensity for gain of a category is less than uniform transitin intensity it labels as "Avoiding".
    transition.behavior.for.gain[[i]] <- matrix(nrow = nrow(transition.intensity[[i]]), ncol = ncol(transition.intensity[[i]]))
    for(j in 1:ncol(transition.intensity[[i]])){
      transition.behavior.for.gain[[i]][,j] <- t(ifelse(abs(t(transition.intensity[[i]][,j]) - uniform.transition[[i]][j]) < Epsilon, "Uniform", ifelse (t(transition.intensity[[i]][,j]) < uniform.transition[[i]][j], "Avoid", "Target")))
    }
    # Make the main diagonal elements of the matrix as "NA".
    # assign row and column names of the matrix.
    diag(transition.behavior.for.gain[[i]]) <- NA
    rownames(transition.behavior.for.gain[[i]]) <- paste(categories)
    colnames(transition.behavior.for.gain[[i]]) <- paste("Transition Behavior for Gain of Category ", categories)

    # Find the annual size of transition from each category to the gaining category.
    Annual.transition.size[[i]] <- matrix(nrow = nrow(transition.intensity[[i]]), ncol = ncol(transition.intensity[[i]]))
    for(z in 1:ncol(transition.intensity[[i]])){
      Annual.transition.size[[i]][,z] <- crosstabulation[[i]][ ,z] / parameters$duration[i]
    }

    diag(Annual.transition.size[[i]]) <- NA
    rownames(Annual.transition.size[[i]]) <- paste(categories)
    colnames(Annual.transition.size[[i]]) <- paste("Transition Size for Gain of Category ", categories)

    # ggplot for transition level
    for (n in 1: length(categories)){
      edited.categories <- categories[-n]
      Annual.transition.size.edited <- Annual.transition.size[[i]][-n, n]
      out.df <- as.data.frame(cbind(edited.categories, Annual.transition.size.edited))
      colnames(out.df) <- c("Categories","transition size")

      edited.transition.intensity <- as.vector(na.omit(transition.intensity[[i]][ ,n]))
      edited.transition.behavior.for.gain <- as.vector(na.omit(transition.behavior.for.gain[[i]][ ,n]))
      edited.uniform.transition <- as.vector(rep(uniform.transition[[i]][n],each= length(categories)-1))
      out.df2 <- as.data.frame(cbind(edited.categories, edited.transition.intensity, edited.uniform.transition, edited.transition.behavior.for.gain))
      out.df2$edited.transition.intensity <- as.numeric(as.character(out.df2$edited.transition.intensity))
      out.df2$edited.uniform.transition <- as.numeric(as.character(out.df2$edited.uniform.transition))
      colnames(out.df2) <- c("Categories", "Transition intensity", "Uniform transition intensity")

      transition.interval <- paste(as.character(parameters$initial.times[i]),
                                   as.character("-"),as.character(parameters$subsequent.times[i]),sep =" ")

      plot(ggplot(out.df, aes(x = edited.categories, y = Annual.transition.size.edited, fill = edited.categories)) +
             geom_bar(width = 0.5, position = position_dodge(width=0.6), stat="identity") + theme_bw() + coord_flip() +
             labs(title = paste("Annual Transition Size for Gain of ", categories[n], "during", transition.interval), x= "Losing Category", y= "Annual Transition Size (# of elements)" ) +
             scale_fill_manual(values = rep("#010158", length(edited.categories))) +
             theme(plot.title = element_text(family = "Times", color="#353535", face="bold", size=14, hjust=0.5)) +
             theme(legend.position = "none", legend.title = element_blank()))

      plot(ggplot(out.df2, aes(x = edited.categories, y = edited.transition.intensity * 100, fill = edited.categories)) +
             geom_bar(width = 0.5, position = position_dodge(width=0.6), stat="identity") + theme_bw() + coord_flip() +
             labs(title = paste("Annual Transition Intensity for Gain of ", categories[n], "during", transition.interval), x = "Losing Category", y = "Annual Transition Intensity (% of category at initial time)") +
             scale_fill_manual(values = rep("#4db8ff", length(edited.categories))) +
             geom_hline(aes(yintercept = out.df2[1,3] * 100), linetype= "dashed", colour = "red", size = 1, show.legend = FALSE) + coord_flip() +
             geom_text(aes(x = length(edited.categories), y = edited.uniform.transition[1] * 100, label = "Target"), family = "Times", colour = "black", angle = 90, vjust = 1.2, show.legend = NA) +
             geom_text(aes(x = length(edited.categories), y = edited.uniform.transition[1] * 100, label = "Avoid"), 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 = "none", legend.title = element_blank()))
    }
  }

  return(list(Annual.Transition.Size = Annual.transition.size, Annual.Transition.Intensity = transition.intensity, Uniform.Transition = uniform.transition, Transition.Behavior.for.Gain = transition.behavior.for.gain))
}