-
Notifications
You must be signed in to change notification settings - Fork 10
/
linearTrend.R
130 lines (113 loc) · 4.19 KB
/
linearTrend.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
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
#' Linear Trend in Remote Sensing
#'
#' Linear trend is useful for mapping forest degradation, land degradation, etc.
#' This algorithm is capable of obtaining the slope of an ordinary least-squares
#' linear regression and its reliability (p-value).
#'
#' @details
#' Linear regression is widely used to analyze forest degradation or land degradation.
#' Specifically, the slope and its reliability are used as main parameters and they
#' can be obtained with this function. On the other hand, logistic regression allows
#' obtaining a degradation risk map, in other words, it is a probability map. Please see
#' the references or see [Tarazona and Maria-Miyasiro (2020)](https://doi.org/10.1016/j.rsase.2020.100337).
#'
#' @references
#' Tarazona, Y., Maria, Miyasiro-Lopez. (2020). Monitoring tropical forest degradation using
#' remote sensing. Challenges and opportunities in the Madre de Dios region, Peru. Remote
#' Sensing Applications: Society and Environment, 19, 100337.
#'
#' Wilkinson, G.N., Rogers, C.E., 1973. Symbolic descriptions of factorial models for
#' analysis of variance. Appl. Stat. 22, 392-399.
#'
#' Chambers, J.M., 1992. Statistical Models in S. CRS Press.
#'
#' @importFrom dplyr bind_rows
#'
#' @param x RasterStack, RasterBrick or Matrix (row*col, n_images)
#' @param type There are two options: "lm" for a linear regression and "glm" for Generalized Linear Models
#' @param ... For "lm" and "glm": arguments to be used to form the default control argument
#' if it is not supplied directly. See \link[stats]{lm} and \link[stats]{glm}.
#'
#' @export
#'
#' @examples
#' library(ForesToolboxRS)
#' library(raster)
#'
#' data(serie_pv)
#'
#' e <- extent(350420.9, 352028.8, -1417869, -1416288)
#' imgs <- crop(serie_pv, e)
#'
#' trend <- linearTrend(x = imgs)
#' plot(trend[[1]]) # raster of slope
#' plot(trend[[2]]) # raster of p-value
#'
linearTrend <- function(x, type = "lm", ...) {
if (type == "lm") {
if (inherits(x, "RasterStack") | inherits(x, "RasterBrick")) {
mat <- raster::as.matrix(x)
} else if (inherits(x, "matrix")) {
mat <- x
} else {
stop(class(x), " class is not supported", call. = TRUE)
}
x_axis <- 1:dim(mat)[2]
my_lms <- lapply(1:dim(mat)[1], function(i) reg <- stats::lm(mat[i, ] ~ x_axis, ...))
summy <- suppressWarnings(
lapply(my_lms, function(x) summary(x)$coefficients[2, 1:4])
)
summy <- as.data.frame(bind_rows(summy))
if (inherits(x, "RasterStack") | inherits(x, "RasterBrick")) {
slope <- raster(x[[1]])
p_value <- raster(x[[2]])
slope[] <- summy[, 1]
p_value[] <- summy[, 4]
result <- raster::stack(slope, p_value)
names(result) <- c("slope", "p_value")
return(result)
} else if (inherits(x, "matrix")) {
result <- array(c(summy[, 1], summy[, 4]), c(dim(x)[1], dim(x)[2], 2))
return(result)
}
} else if (type == "glm") {
if (inherits(x, "RasterStack") | inherits(x, "RasterBrick")) {
names(x[[1]]) <- c("FD")
df <- raster::as.matrix(x)
mod <- stats::glm(FD ~ ., data = df, ...)
prob <- predict(mod, type = "response")
result <- raster(x[[1]])
result[] <- prob
names(result) <- c("MapProbability")
return(result)
} else if (inherits(x, "list")) {
if (class(try(raster::stack(x), silent = TRUE)) == "try-error") {
master <- raster(x[1])
meta <- raster(extent(master), nrows = master@nrows, ncols = master@ncols)
varib_rast <- c()
for (k in 1:length(x)) {
ras_ini <- raster(x[[k]])
ras_post <- raster::resample(ras_ini, meta, method = "ngb")
ras_post[is.na(ras_post)] <- 0
varib_rast <- c(varib_rast, ras_post)
}
x <- varib_rast
names(x[[1]]) <- c("FD")
} else {
x <- raster::stack(x)
names(x[[1]]) <- c("FD")
}
df <- raster::as.matrix(x)
mod <- stats::glm(FD ~ ., data = df, ...)
prob <- predict(mod, type = "response")
result <- master
result[] <- prob
names(result) <- c("MapProbability")
return(result)
} else {
stop(class(x), " class is not supported", call. = TRUE)
}
} else {
stop(" Method is not supported", call. = TRUE)
}
}