R/chg_traj.R

In azvoleff/teamlucc: TEAM land use and cover change data processing toolkit

#' Calculate Quantum GIS format color map file
#'
#' This function will format a lookup table (lut) into a QGIS color map 
#' suitable for import into Quantum GIS (QGIS) for color coding a raster file 
#' using the "singleband pseudocolor" band rendering option.
#'
#' @export
#' @import raster
#' @param x a lookup table as output by \code{traj_lut}
#' @param out_file filename for the output file
#' @param color_codes an (optional) list of colors to use in the color table. Should 
#' have the same length as the number of rows in x.
#' @param minlength number of characters used to abbreviate t0_name and t1_name 
#' columns when building transition names (only used it t0_name and t1_name are 
#' present in \code{x})
#' @examples
#' lut <- traj_lut(c(1, 2), c("NonForest", "Forest"))
qgis_colormap <- function(x, out_file, color_codes, minlength=5) {
    if (extension(out_file) == "txt") {
        stop("out_file must be a text file")
    }
    if (file_test('-f', out_file)) {
        stop("out_file already exists")
    }

    if (missing(color_codes)) {
        gg_color_hue <- function(n) {
            # From: http://bit.ly/1yXZOlv
            hues = seq(15, 375, length=n+1)
            hcl(h=hues, l=65, c=100)[1:n]
        }
        color_codes <- gg_color_hue(nrow(x))
    }

    color_codes <- t(col2rgb(color_codes))
    x$r <- color_codes[, 1]
    x$g <- color_codes[, 2]
    x$b <- color_codes[, 3]

    if ("t0_name" %in% names(x) & ("t1_name" %in% names(x))) {
        x$trans <- apply(cbind(abbreviate(x$t0_name, minlength),
                               abbreviate(x$t1_name, minlength)),
                         1, paste, collapse="->")
    } else {
        x$trans <- apply(cbind(x$t0_code, x$t1_code), 1, paste, collapse="->")
    }

    cat(c("# Color map generated by teamlucc", "INTERPOLATION:INTERPOLATED"), sep="\n", file=out_file)

    out_data <- apply(cbind(x$Code, x$r, x$g, x$b, 255, x$trans), 1, paste, 
                      collapse=",")

    cat(out_data, sep="\n", file=out_file, append=TRUE)
}

#' Calculate change-trajectory lookup table
#'
#' This function will format a lookup table (lut) to allow coding change 
#' trajectories. Useful for use in conjunction with \code{\link{chg_traj}}.
#'
#' @export
#' @param class_codes a list of integer codes used to code land use/cover 
#' classes
#' @param class_names an (optional) list of class names as character vectors
#' @examples
#' lut <- traj_lut(c(1, 2), c("NonForest", "Forest"))
traj_lut <- function(class_codes, class_names=NULL) {
    lut <- expand.grid(t0_code=class_codes, t1_code=class_codes)
    if (!is.null(class_names)) {
        if (length(class_names) != length(class_codes)) {
            stop('class_names must be NULL or a vector of length equal to number of classes in initial image')
        }
        lut$t0_name <- class_names[match(lut$t0_code, class_codes)]
        lut$t1_name <- class_names[match(lut$t1_code, class_codes)]
    }
    # Code trajectories by summing t0 and t1 after multiplying t1 by the number 
    # of classes.
    lut$Code <- lut$t0_code + lut$t1_code * length(class_codes)
    # Exclude classes that are persistence - CVAPS doesn't directly code the 
    # class for classes that persist
    lut <- lut[!(lut$t0_code == lut$t1_code), ]
    return(lut)
}

#' Calculate change-trajectory image
#'
#' This function will calculate trajectories of land cover change using the 
#' Change Vector Analysis in Posterior Probability Space (CVAPS) approach of 
#' comparing posterior probabilities of class membership with an automatically 
#' determined threshold. Areas of no change are coded as -1. A lookup table for 
#' the codes output by \code{chg_traj} can be calculated with \code{traj_lut}.
#'
#' This function will run in parallel if a parallel backend is registered with 
#' \code{\link{foreach}}.
#'
#' @export
#' @importFrom spatial.tools rasterEngine
#' @param chg_mag change magnitude \code{RasterLayer} from \code{CVAPS}
#' @param chg_dir change direction \code{RasterLayer} from \code{CVAPS}
#' @param chg_threshold the threshold to use as a minimum when determining change 
#' areas (can use \code{DFPS} to determine this value).
#' @param filename filename to save the output \code{RasterLayer} to disk 
#' (optional)
#' @param overwrite whether to overwrite existing files (otherwise an error 
#' will be raised)
#' @param ... additional parameters to pass to rasterEngine
#' @return a {RasterLayer} of change trajectories, with change trajectories 
#' coded as in the \code{lut} output by \code{traj_lut}
#' @references Chen, J., P. Gong, C.  He, R.  Pu, and P.  Shi.  2003.
#' Land-use/land-cover change detection using improved change-vector analysis.
#' Photogrammetric Engineering and Remote Sensing 69:369-380.
#' 
#' Chen, J., X. Chen, X. Cui, and J. Chen. 2011. Change vector analysis in
#' posterior probability space: a new method for land cover change detection.
#' IEEE Geoscience and Remote Sensing Letters 8:317-321.
#' @examples
#' \dontrun{
#' t0_train_data <- get_pixels(L5TSR_1986, L5TSR_1986_2001_training, "class_1986",training=.6)
#' t0_model <- train_classifier(t0_train_data)
#' t0_preds <- classify(L5TSR_1986, t0_model)
#' t1_train_data <- get_pixels(L5TSR_2001, L5TSR_1986_2001_training, "class_2001", training=.6)
#' t1_model <- train_classifier(t1_train_data)
#' t1_preds <- classify(L5TSR_2001, t1_model)
#' t0_t1_chgmag <- chg_mag(t0_preds$probs, t1_preds$probs)
#' t0_t1_chgdir <- chg_dir(t0_preds$probs, t1_preds$probs)
#' 
#' lut <- traj_lut(t0_preds$codes$code, t0_preds$codes$class)
#' t0_t1_chgtraj <- chg_traj(lut, t0_t1_chgmag, t0_t1_chgdir, .5)
#' 
#' # Change areas are coded following the above lookup-table (lut):
#' plot(t0_t1_chgtraj)
#' 
#' # No change areas are -1:
#' plot(t0_t1_chgtraj == -1)
#' }
chg_traj <- function(chg_mag, chg_dir, chg_threshold, filename, 
                     overwrite=FALSE, ...) {
    if (nlayers(chg_mag) > 1) stop('chg_mag has more than 1 layer')
    if (nlayers(chg_dir) > 1) stop('chg_dir has more than 1 layer')
    compareRaster(chg_mag, chg_dir)
    if (!missing(filename) && file_test('-f', filename) && !overwrite) {
        stop('output file already exists and overwrite=FALSE')
    }

    calc_chg_traj <- function(chg_mag, chg_dir, chg_threshold, ...) {
        # Trajectories in chg_dir were coded by summing t0 and t1 classes after 
        # multiplying t1 class by the number of classes
        chg_dir[chg_mag < chg_threshold] <- -1
        chg_dir[is.na(chg_dir)] <- -2
        chg_dir <- array(chg_dir, dim=c(dim(chg_mag)[1], dim(chg_mag)[2], 1))
        return(chg_dir)
    }
    out <- rasterEngine(chg_mag=chg_mag, chg_dir=chg_dir, fun=calc_chg_traj,
                        args=list(chg_threshold=chg_threshold),
                        datatype='INT2S', ...)

    # spatial.tools doesn't properly handle NA values for integer layers, so 
    # they were coded as -2 above - now recode them as NA
    out[out == -2] <- NA

    # spatial.tools can only output the raster package grid format - so output 
    # to a tempfile in that format then copy over to the requested final output 
    # format if a filename was supplied
    if (!missing(filename)) {
        out <- writeRaster(out, filename=filename, overwrite=overwrite, 
                           datatype='INT2S')
    }

    return(out)
}