Nothing
# This file is a part of zonator package
# Copyright (C) 2012-2014 Joona Lehtomaki <joona.lehtomaki@gmai.com>. All rights
# reserved.
# This program is open source software; you can redistribute it and/or modify
# it under the terms of the FreeBSD License (keep this notice):
# http://en.wikipedia.org/wiki/BSD_licenses
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
#' Compare matrices in various ways.
#'
#' Function can be used to compare two Zonation output rasters with one of the
#' following functions (part of zonator package):
#' \itemize{
#' \item{correlation}{}
#' \item{substraction}{}
#' \item{ferquency}{(NOT IMPLEMENTED)}
#' \item{coverage}{}
#' }
#'
#' @keywords post-processing ppa
#' @author Joona Lehtomaki <joona.lehtomaki@@gmail.com>
#'
#' @param x Numeric matrix.
#' @param y Numeric matrix.
#' @param fun Function used for the numeric comparison.
#' @param ... Further arguments passed on to selected comparison function.
#'
#' @return A DataFrame with each row containg columns title, count, and catid
#'
#' @export
#' @seealso correlation substraction frequency coverage
comp <- function(x, y, fun="correlation", ...) {
# Check the data
if (!is.matrix(x) | !is.matrix(y)) {
stop("Both inputs must be matrices.")
}
if (dim(x)[1] != dim(y)[1] | dim(x)[2] != dim(y)[2]) {
stop("Matrix dimensions do not match.")
}
switch(fun,
correlation = correlation(x, y, ...),
substraction = (x - y),
# NOT IMPLEMENTED
# frequency = selection.frequency(x, y, ...),
coverage = selection_coverage(x, y, ...))
}
#' Correlation between two matrices.
#'
#' Calculate correlation between two matrices using \code{\link{cor}}. A
#' group of specific threshold can be set, in which case the correlations are
#' calculated incrementally for values above the thresholds.
#'
#' @keywords post-processing ppa
#' @author Joona Lehtomaki <joona.lehtomaki@@gmail.com>
#'
#' @param x Numeric matrix.
#' @param y Numeric matrix.
#' @param method Character string correlation method used (default: 'kendall').
#' @param thresholds Numeric vector of thresholds used (default: c(0)).
#'
#' @return A list with 2 items:
#' \item{thresholds}{Correlations between 2 matrices with values above a given threshold.}
#' \item{total}{Overall correlation between the 2 matrices.}
#'
#' @importFrom stats cor
#'
#' @export
#' @seealso \code{\link{cor}}
correlation <- function(x, y, method='kendall', thresholds=c(0)) {
res <- c()
for (i in 1:length(thresholds)) {
x.sel <- as.vector(x[which(x > thresholds[i])])
y.sel <- as.vector(y[which(y > thresholds[i])])
res <- append(res, cor(x.sel, y.sel, method=method))
}
res <- data.frame(res, row.names=thresholds)
colnames(res) <- "correlation"
return(list(thresholds=res, total=cor(as.vector(x), as.vector(y),
method=method)))
}
#' Intersection of two coverages.
#'
#' Calculate how much two coverages (as defined by values greater than a given
#' threshold in two numeric matrices) overlap.
#'
#' @keywords post-processing ppa
#' @author Joona Lehtomaki <joona.lehtomaki@@gmail.com>
#'
#' @param x Numeric matrix.
#' @param y Numeric matrix.
#' @param thresholds Numeric vector of thresholds used.
#'
#' @return A list with 2 items:
#' \item{thresholds}{Correlations between 2 matrices with values above a given threshold.}
#' \item{total}{Overall correlation between the 2 matrices.}
#'
#' @export
selection_coverage <- function(x, y, thresholds) {
covs <- c()
total <- c()
for (thresh in thresholds) {
sel1 <- which(x >= thresh)
sel2 <- which(y >= thresh)
# What fraction of the original coverage x in selection?
total <- append(total, length(sel1) / length(x))
covs <- append(covs, sum(sel1 %in% sel2) / length(sel1))
}
res <- data.frame(total=total, cover=covs, row.names=thresholds)
return(res)
}
#' Calculate the Jaccard coefficient.
#'
#' The Jaccard coefficient measures similarity between sample sets, and is
#' defined as the size of the intersection divided by the size of the union of
#' the sample sets. The Jaccard coefficient can be calculated for a subset of
#' rasters provided by using the threshold argument.
#'
#' Min and max values must be provided for both RasterLayer objects \code{x}
#' and \code{y}. Method can be used with RasterLayers of any value range, but
#' the defaults [0.0, 1.0] are geared towards comparing Zonation rank priority
#' rasters. Limits provided are inclusive.
#'
#' @keywords post-processing ppa
#' @author Joona Lehtomaki <joona.lehtomaki@@gmail.com>
#'
#' @param x raster object.
#' @param y raster object.
#' @param x.min Numeric minimum threshold value for \code{x} to be used
#' (default 0.0).
#' @param x.max Numeric maximum threshold value for \code{x} to be used
#' (default 1.0).
#' @param y.min Numeric minimum threshold value for \code{y} to be used
#' (default 0.0).
#' @param y.max Numeric maximum threshold value for \code{y} to be used
#' (default 1.0).
#' @param warn.uneven Logical indicating whether a warning is raised if the
#' compared raster coverages are very (>20x) uneven.
#' @param limit.tolerance integer values that defines to which precision x and
#' y limits are rounded to. This helps e.g. with values that close to 0 but
#' not quite 0 (default: 4, i.e. round(x, 4)).
#' @param disable.checks logical indicating if the input limit values are
#' checked against the actual raster values in \code{x} and \code{y}.
#'
#' @return A numeric value [0, 1]
#'
#' @export
#'
jaccard <- function(x, y, x.min=0.0, x.max=1.0, y.min=0.0, y.max=1.0,
warn.uneven=FALSE, limit.tolerance=4,
disable.checks=FALSE) {
if (!disable.checks) {
# Check the input values
x.min.value <- round(cellStats(x, stat="min"), limit.tolerance)
x.max.value <- round(cellStats(x, stat="max"), limit.tolerance)
y.min.value <- round(cellStats(y, stat="min"), limit.tolerance)
y.max.value <- round(cellStats(y, stat="max"), limit.tolerance)
if (x.min < x.min.value) {
stop(paste0("Minimum threshold value for x ("), x.min, ") smaller than
the computed minimum value in x (", x.min.value, ")")
}
if (x.max > x.max.value) {
stop(paste0("Maximum threshold value for x ("), x.max, ") smaller than
the computed maximum value in x (", x.max.value, ")")
}
if (x.min >= x.max) {
stop(paste0("Minimum threshold value for x ("), x.min, ") smaller than
maximum threshold value for x (", x.max, ")")
}
if (y.min < y.min.value) {
stop(paste0("Minimum threshold value for y ("), y.min, ") smaller than
the computed minimum value in y (", y.min.value, ")")
}
if (y.max > y.max.value) {
stop(paste0("Maximum threshold value for y ("), y.max, ") smaller than
the computed maximum value in y (", y.max.value, ")")
}
if (y.min >= y.max) {
stop(paste0("Minimum threshold value for y ("), y.min, ") smaller than
maximum threshold value for y (", y.max, ")")
}
# Comparisons using just the defaults is probably not feasible
if (x.min == 0.0 & x.max == 1.0 & y.min == 0.0 & y.max == 1.0) {
warning("Using all the defaults for x and y ranges")
}
} else {
message("Input limit checks disabled")
}
# [fixme] - using cellStats(X, "sum") should be safe as we're dealing with
# binary 0/1 rasters. count() would be preferable, but apparently raster
# (>= 2.2 at least) doesn't support it anymore.
# Get the values according to the limits provided
x.bin <- (x >= x.min & x <=x.max)
y.bin <- (y >= y.min & y <=y.max)
if (warn.uneven) {
x.size <- cellStats(x.bin, "sum")
y.size <- cellStats(y.bin, "sum")
# Sort from smaller to larger
sizes <- sort(c(x.size, y.size))
if (sizes[2] / sizes[1] > 20) {
warning("The extents of raster values above the threshhold differ more",
"than 20-fold: Jaccard coefficient may not be informative.")
}
}
# Calculate the intersection of the two rasters, this is given by adding
# the binary rasters together -> 2 indicates intersection
combination <- x.bin + y.bin
intersection <- combination == 2
# Union is all the area covered by the both rasters
union <- combination >= 1
return(cellStats(intersection, "sum") / cellStats(union, "sum"))
}
#' Calculate Jaccard coefficients bewteen all the RasterLayers within a single
#' RasterStack.
#'
#' This method is a utility method that is intented to be used to compare
#' top-fractions of the landscape. Thus, x.max and y.max for
#' \code{\link{jaccard}} are fixed to 1.0.
#'
#' @keywords post-processing ppa
#' @author Joona Lehtomaki <joona.lehtomaki@@gmail.com>
#'
#' @param stack RasterStack-object.
#' @param thresholds Numeric vector values of thresholds.
#' @param ... additional arguments passed on to \code{\link{jaccard}}.
#'
#' @return Dataframe with Jaccard coefficients between all the RasterLayers.
#'
#' @export
#' @seealso \code{\link{jaccard}}
cross_jaccard <- function(stack, thresholds, ...) {
all.jaccards <- list()
for (threshold in thresholds) {
jaccards <- matrix(nrow=nlayers(stack), ncol=nlayers(stack))
for (i in 1:nrow(jaccards)) {
for (j in 1:ncol(jaccards)) {
if (i == j) {
jaccards[i, j] <- 1
} else {
# See the complement, if it's not NA then the pair has already been
# compared
if (is.na(jaccards[j, i])) {
message(paste0("Calculating Jaccard index for [", threshold, ", ",
1.0, "] between ", names(stack[[i]]), " and ",
names(stack[[j]])))
jaccards[i, j] <- jaccard(stack[[i]], stack[[j]],
x.min=threshold, x.max=1.0,
y.min=threshold, y.max=1.0, ...)
} else {
jaccards[i, j] <- jaccards[j, i]
}
}
}
}
jaccards <- as.data.frame(jaccards)
colnames(jaccards) <- names(stack)
rownames(jaccards) <- names(stack)
all.jaccards[[as.character(threshold)]] <- jaccards
}
return(all.jaccards)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.