-
Notifications
You must be signed in to change notification settings - Fork 10
/
fusionRS.R
133 lines (118 loc) · 4.46 KB
/
fusionRS.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
131
132
133
#' Fusion of images with different observation geometries
#'
#' This algorithm allows to fusion images coming from different spectral sensors
#' (e.g., optical-optical, optical and SAR, or SAR-SAR). It is also possible to obtain the contribution (%) of each variable
#' in the fused image.
#'
#' @return
#' A list with Fused images (RasterStack), Variance, Proportion of variance, Cumulative variance,
#' Correlation and Contribution in percentage.
#'
#' @note
#' Before executing the function, it is recommended that images coming from different
#' sensors or from the same sensor have a co-registration.
#'
#' @importFrom stats prcomp na.omit
#' @importFrom raster getValues as.data.frame brick extent plotRGB
#' @importFrom factoextra get_pca_var
#'
#' @param x Optical image. It should be RasterStack or RasterBrick
#' @param y Radar image. It should be RasterStack or RasterBrick
#' @param stand.varb Logical. If \code{TRUE}, the PCA is calculated using the correlation
#' matrix (standardized variables) instead of the covariance matrix (non-standardized variables)
#' @param na.rm If \code{TRUE} the NA values of the images will be omitted from the analysis
#' @param verbose This parameter is Logical. It prints progress messages during execution
#'
#' @examples
#' \dontrun{
#' library(ForesToolboxRS)
#' library(raster)
#'
#' # Load example datasets
#' data(img_optical)
#' data(img_radar)
#'
#' # Fusion of images - optical and radar data
#' fusion <- fusionRS(x = img_optical, y = img_radar, na.rm = TRUE)
#' plotRGB(fusion[[1]], 1, 2, 3, axes = FALSE, stretch = "lin", main = "Fused images")
#' }
#' @export
#'
fusionRS <- function(x, y, stand.varb = TRUE, na.rm = FALSE, verbose = FALSE) {
if (verbose) {
message(paste0(paste0(rep("*", 10), collapse = ""), " Verifying the same extent ", paste0(rep("*", 10), collapse = "")))
}
if (inherits(x, "RasterStack") | inherits(x, "RasterBrick") & inherits(y, "RasterStack") | inherits(y, "RasterBrick")) {
if (extent(x) == extent(y)) {
img <- raster::stack(x, y)
} else {
stop(" The extent of the images are different.", call. = TRUE)
}
} else {
stop(c(class(x), class(y)), " This classes are not supported yet.", call. = TRUE)
}
df <- as.data.frame(img)
# Omit NA
if (na.rm) {
df <- na.omit(df)
}
# standardized variables?
if (stand.varb) {
if (verbose) {
message(paste0(paste0(rep("*", 10), collapse = ""), " Standardizing variables ", paste0(rep("*", 10), collapse = "")))
}
acp <- prcomp(df, center = TRUE, scale = TRUE) # standardized variables
} else {
acp <- prcomp(df, center = FALSE, scale = FALSE) # non-standardized variables
}
if (verbose) {
message(paste0(paste0(rep("*", 10), collapse = ""), " Contributions and correlations ", paste0(rep("*", 10), collapse = "")))
}
var <- summary(acp)$importance[1, ]^2 # Variance
pov <- summary(acp)$importance[2, ] # Proportion of variance
varAcum <- summary(acp)$importance[3, ] # Cumulative variance
corr <- get_pca_var(acp)$cor # Correlation
contri <- get_pca_var(acp)$contrib # Contribution in %
for (i in 1:dim(df)[2]) {
colnames(corr)[i] <- paste("PC", i, sep = "")
rownames(corr)[i] <- paste("Band", i, sep = "")
colnames(contri)[i] <- paste("PC", i, sep = "")
rownames(contri)[i] <- paste("Band", i, sep = "")
}
# Storing in a raster the principal components
acpY <- img
vr <- getValues(acpY)
p <- which(!is.na(vr)) # Positions of NA
# We assign a raster format to each PC
for (k in 1:dim(df)[2]) {
acpY[[k]][p] <- acp$x[, k]
}
results <- list(FusionRS = acpY, var = var, pov = pov, varAcum = varAcum, corr = corr, contri = contri)
names(results) <- c(
"Fused_images", "Variance", "Proportion_of_variance",
"Cumulative_variance", "Correlation", "Contribution_in_pct"
)
return(structure(results, class = "fusionRS"))
}
#' Print of the "fusionRS" class
#'
#' @param x Object of class "fusionRS".
#' @param ... Not used.
#'
#' @export
#'
print.fusionRS <- function(x, ...) {
cat("******************** ForesToolboxRS - FUSION OF IMAGES ********************\n")
cat("\n**** Fused images ****\n")
methods::show(x$Fused_images)
cat("\n**** Variance ****\n")
print(x$Variance)
cat("\n**** Proportion_of_variance ****\n")
print(x$Proportion_of_variance)
cat("\n**** Cumulative_variance ****\n")
print(x$Cumulative_variance)
cat("\n**** Correlation ****\n")
print(x$Correlation)
cat("\n**** Contribution_in_percentage ****\n")
print(x$Contribution_in_pct)
}