-
Notifications
You must be signed in to change notification settings - Fork 29
/
rasterCorrelation.R
66 lines (66 loc) · 2.34 KB
/
rasterCorrelation.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
#' @title Raster correlation
#' @description
#' Performs a moving window correlation between two rasters
#'
#' @param x A terra SpatRaster class object for x
#' @param y A terra SpatRasterclass object for y
#' @param s Scale of window. Can be a single value, two
#' values for uneven window or a custom matrix.
#' Must be odd number (eg., s=3, for 3x3 window or
#' s=c(3,5) for 3 x 5 window)
#' @param type Type of output, options are: "pearson", "spearman",
#' "covariance"
#'
#' @note
#' The NA behavior is set to na.rm = TRUE to make default outputs
#' consistent between the terra and raster packages.
#'
#' @return A terra SpatRaster class object
#'
#' @author Jeffrey S. Evans <jeffrey_evans@@tnc.org>
#'
#' @examples
#' \donttest{
#' library(terra)
#'
#' r <- rast(system.file("ex/logo.tif", package="terra"))
#' x <- r[[1]]
#' y <- r[[3]]
#'
#' r.cor <- rasterCorrelation(x, y, s = 5, type = "spearman")
#' plot(r.cor)
#' }
#'
#' @export rasterCorrelation
rasterCorrelation <- function(x, y, s = 3, type = "pearson") {
if (!inherits(x, "SpatRaster"))
stop(deparse(substitute(x)), " must be a terra SpatRaster object")
if (!inherits(y, "SpatRaster"))
stop(deparse(substitute(y)), " must be a terra SpatRaster object")
if(!inherits(s, "matrix")) {
if(length(s) > 2)
stop( "Specified window exceeds 2 dimensions")
if(length(s) == 1) { s = rep(s, 2) }
s <- matrix(1,s[1],s[2])
}
tmpXY = x * y
xBar <- terra::focal(x, w = s, fun = mean)
yBar <- terra::focal(y, w = s, fun = mean)
xyBar <- terra::focal(tmpXY, w = s, fun = mean)
if(type == "covariance") {
return( xyBar - (xBar * yBar) )
} else if (type == "pearson") {
coVar <- xyBar - (xBar * yBar)
xStd <- terra::focal(x, w = s, fun=stats::sd)
yStd <- terra::focal(y, w = s, fun=stats::sd)
xyStd <- xStd * yStd
return( coVar / xyStd )
} else if (type == "spearman") {
coVar <- xyBar - (xBar * yBar)
xStd <- terra::focal(x, w = s, fun=stats::sd)
yStd <- terra::focal(y, w = s, fun=stats::sd)
return( coVar / (xStd * yStd) )
} else {
stop( "Only pearson, spearman correlation or covariance are supported")
}
}