Skip to content

Commit e0eff23

Browse files
jbferetjbferet
committed
# biodivMapR2 v2.3.11
## change - use continuumRemoval from package prospectr instead of internal function # biodivMapR2 v2.3.10 ## change - add option "BIGTIFF=IF_SAFER" when calling sf::gdal_utils to produce mosaics
1 parent f825637 commit e0eff23

File tree

10 files changed

+287
-84
lines changed

10 files changed

+287
-84
lines changed

‎DESCRIPTION‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
Package: biodivMapR
22
Title: biodivMapR: an R package for a- and ß-diversity mapping using remotely-sensed images
3-
Version: 2.3.10
3+
Version: 2.3.11
44
Authors@R: c(person(given = "Jean-Baptiste",
55
family = "Feret",
66
email = "jb.feret@teledetection.fr",
@@ -44,7 +44,7 @@ Imports:
4444
tools,
4545
tidyr,
4646
vegan,
47-
crsuggest
47+
crsuggest,
4848
Suggests:
4949
ggplot2,
5050
gridExtra,

‎NAMESPACE‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -19,7 +19,9 @@ export(compute_bc_diss)
1919
export(compute_mask_iqr)
2020
export(compute_mask_iqr_tiles)
2121
export(compute_nn_from_ordination)
22+
export(continuumRemoval)
2223
export(continuum_removal)
24+
export(crfun)
2325
export(define_pixels_per_iter)
2426
export(exclude_spectral_domains)
2527
export(explore_kmeans)
@@ -125,6 +127,7 @@ importFrom(future,plan)
125127
importFrom(future,sequential)
126128
importFrom(future.apply,future_lapply)
127129
importFrom(future.apply,future_mapply)
130+
importFrom(grDevices,chull)
128131
importFrom(methods,as)
129132
importFrom(parallel,makeCluster)
130133
importFrom(parallel,stopCluster)
@@ -144,6 +147,7 @@ importFrom(sf,st_sf)
144147
importFrom(sf,st_sfc)
145148
importFrom(snow,splitRows)
146149
importFrom(stats,IQR)
150+
importFrom(stats,approx)
147151
importFrom(stats,as.dist)
148152
importFrom(stats,cmdscale)
149153
importFrom(stats,cor.test)
@@ -154,6 +158,7 @@ importFrom(stats,prcomp)
154158
importFrom(stats,quantile)
155159
importFrom(stats,runif)
156160
importFrom(stats,sd)
161+
importFrom(stats,splinefun)
157162
importFrom(terra,aggregate)
158163
importFrom(terra,blocks)
159164
importFrom(terra,buffer)

‎NEWS.md‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,7 @@
1+
# biodivMapR2 v2.3.11
2+
## change
3+
- use continuumRemoval from package prospectr instead of internal function
4+
15
# biodivMapR2 v2.3.10
26
## change
37
- add option "BIGTIFF=IF_SAFER" when calling sf::gdal_utils to produce mosaics

‎R/apply_continuum_removal.R‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,16 @@ apply_continuum_removal <- function(spectral_data, spectral) {
1818
# avoids memory crash
1919
max_nb_values <- 2e6
2020
nb_cr <- ceiling(nb_values / max_nb_values)
21+
22+
# # ensure first wavelength is not 0
23+
# if (max(spectral_data[,1])>1)
24+
# offset <- 1
25+
# if (max(spectral_data[,1])<1)
26+
# offset <- 1e-4
27+
# sel <- which(spectral_data[,1]==0)
28+
# if (length(sel)>0)
29+
# spectral_data[sel,1] <- spectral_data[sel,1] + offset
30+
2131
spectral_data <- snow::splitRows(spectral_data, nb_cr)
2232
spectral_data_tmp <- lapply(spectral_data, FUN = continuum_removal,
2333
spectral_bands = spectral$Wavelength)

‎R/continuumRemoval.R‎

Lines changed: 77 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,77 @@
1+
#' @title Continuum Removal
2+
#' @description adapted from continuumRemoval function in R package propsectr
3+
#' https://github.com/l-ramirez-lopez/prospectr/blob/main/R/continuumRemoval.R
4+
#'
5+
#' Compute the continuum removed values of a data matrix or vector
6+
#' @usage
7+
#' continuumRemoval(X, wav, type = c("R", "A"),
8+
#' interpol = c("linear", "spline"),
9+
#' method = c("division", "substraction"))
10+
#' @param X a numeric matrix or vector to process (optionally a data frame that can
11+
#' be coerced to a numerical matrix).
12+
#' @param wav optional. A numeric vector of band positions.
13+
#' @param type the type of data: 'R' for reflectance (default), 'A' for
14+
#' absorbance.
15+
#' @param interpol the interpolation method between points on the convex hull:
16+
#' 'linear' (default) or 'spline'.
17+
#' @param method normalization method: 'division' (default) or 'subtraction'
18+
#' (see details section).
19+
#' @author Antoine Stevens & \href{https://orcid.org/0000-0002-5369-5120}{Leonardo Ramirez-Lopez}
20+
#' @return a matrix or vector with the filtered spectra.
21+
#' @details
22+
#' The continuum removal technique was introduced by Clark and Roush (1984)
23+
#' as a method to highlight energy absorption features of minerals.
24+
#' It can be viewed as a way to perform albedo normalization.
25+
#' The algorithm find points lying on the convex hull (local maxima or envelope)
26+
#' of a spectrum, connects the points by linear or spline interpolation and
27+
#' normalizes the spectrum by dividing (or subtracting) the input data by the
28+
#' interpolated line.
29+
#' @references
30+
#' Clark, R.N., and Roush, T.L., 1984. Reflectance Spectroscopy: Quantitative
31+
#' Analysis Techniques for Remote Sensing Applications. J. Geophys. Res. 89,
32+
#' 6329-6340.
33+
#' @export
34+
35+
continuumRemoval <- function(X,
36+
wav,
37+
type = c("R", "A"),
38+
interpol = c("linear", "spline"),
39+
method = c("division", "substraction")) {
40+
if (is.data.frame(X))
41+
X <- as.matrix(X)
42+
43+
type <- match.arg(type)
44+
interpol <- match.arg(interpol)
45+
method <- match.arg(method)
46+
47+
if (type == "A")
48+
X <- 1 / X
49+
50+
if (is.matrix(X)) {
51+
if (missing(wav))
52+
wav <- seq_len(ncol(X))
53+
if (length(wav) != ncol(X))
54+
stop("length(wav) should be equal to ncol(X)")
55+
cont <- t(apply(X, 1, function(x) crfun(x, wav, interpol)))
56+
} else {
57+
cont <- crfun(X, wav, interpol)
58+
}
59+
60+
if (method == "division") {
61+
cr <- X / cont
62+
} # like ENVI
63+
else {
64+
cr <- 1 + X - cont
65+
}
66+
67+
if (type == "A")
68+
cr <- 1 / cr - 1
69+
70+
if (is.matrix(X)) {
71+
colnames(cr) <- wav
72+
rownames(cr) <- rownames(X)
73+
} else {
74+
names(cr) <- wav
75+
}
76+
return(cr)
77+
}

‎R/continuum_removal.R‎

Lines changed: 85 additions & 80 deletions
Original file line numberDiff line numberDiff line change
@@ -14,90 +14,95 @@ continuum_removal <- function(mat_init, spectral_bands, p = NULL) {
1414

1515
# Filter and prepare data prior to continuum removal
1616
cr_data <- filter_prior_cr(mat_init, spectral_bands)
17-
mat_init <- cr_data$mat_init
18-
nb_bands <- dim(mat_init)[2]
19-
cr_data$mat_init <- c()
20-
spectral_bands <- cr_data$spectral_bands
17+
nb_bands <- dim(cr_data$mat_init)[2]
2118
nb_samples <- cr_data$nb_samples
19+
# mat_init <- cr_data$mat_init
20+
# cr_data$mat_init <- c()
21+
# spectral_bands <- cr_data$spectral_bands
2222
nb_samples_update <- length(cr_data$samples_to_keep)
23+
2324
# if samples to be considered
2425
if (nb_samples > 0) {
25-
# initialization:
26-
# spectral band corresponding to each element of the data matrix
27-
lambda <- repmat(matrix(spectral_bands, nrow = 1), nb_samples_update, 1)
28-
# prepare matrices used to check evolution of the CR process
29-
# - elements still not processed through continuum removal: initialization to 1
30-
still_need_cr <- matrix(1, nrow = nb_samples_update, ncol = nb_bands)
31-
# - value of the convex hull: initially set to 0
32-
convex_hull <- matrix(0, nrow = nb_samples_update, ncol = nb_bands)
33-
# - reflectance value for latest interception with convex hull:
34-
# initialization to value of the first reflectance measurement
35-
intercept_hull <- repmat(matrix(mat_init[, 1], ncol = 1), 1, nb_bands)
36-
# - spectral band of latest interception
37-
latest_intercept <- repmat(X = matrix(spectral_bands[1], ncol = 1),
38-
m = nb_samples_update, n = nb_bands)
39-
# number of spectral bands found as intercept
40-
nb_intercept <- 0
41-
# continues until arbitrary stopping criterion:
42-
# stops when reach last spectral band (all values before last = 0)
43-
# while (max(still_need_cr[, seq_len(nb_bands - 2)]) == 1 & (nb_intercept <= (nb_bands / 2))) {
44-
while (max(still_need_cr[, seq_len((nb_bands - 2))]) == 1) {
45-
nb_intercept <- nb_intercept + 1
46-
# identify samples still needing continuum removal
47-
sel <- which(still_need_cr[,(nb_bands-2)]==1)
48-
# update variables to process samples needing CR only
49-
nb_samples_update_tmp <- length(sel)
50-
lambda_tmp <- lambda[sel,]
51-
mat_init_tmp <- mat_init[sel,]
52-
latest_intercept_tmp <- latest_intercept[sel,]
53-
still_need_cr_tmp <- still_need_cr[sel,]
54-
convex_hull_tmp <- convex_hull[sel,]
55-
intercept_hull_tmp <- intercept_hull[sel,]
56-
# Mstep give the position of the values to be updated
57-
update_data <- matrix(1, nrow = nb_samples_update_tmp, ncol = nb_bands)
58-
update_data[, nb_bands] <- 0
59-
# initial step: first column set to 0; following steps: all bands below
60-
# max of the convex hull are set to 0
61-
update_data[which((lambda_tmp - latest_intercept_tmp) < 0)] <- 0
62-
# compute slope for each coordinate
63-
slope <- as.matrix((mat_init_tmp - intercept_hull_tmp) / (lambda_tmp - latest_intercept_tmp) * still_need_cr_tmp)
64-
# set current spectral band and previous bands to -9999
65-
if (!length(which(still_need_cr_tmp == 0)) == 0) {
66-
slope[which(still_need_cr_tmp == 0)] <- -9999
67-
}
68-
if (!length(which(is.na(slope))) == 0) {
69-
slope[which(is.na(slope))] <- -9999
70-
}
71-
# get max index for each row and convert into linear index
72-
index_max_slope <- RowToLinear(max.col(slope, ties.method = "last"),
73-
nb_samples_update_tmp, nb_bands)
74-
# !!!! OPTIM: replace repmat with column operation
75-
# update coordinates of latest intercept
76-
latest_intercept_tmp <- repmat(matrix(lambda_tmp[index_max_slope],
77-
ncol = 1), 1, nb_bands)
78-
# update latest intercept
79-
intercept_hull_tmp <- repmat(matrix(as.matrix(mat_init_tmp)[index_max_slope],
80-
ncol = 1), 1, nb_bands)
81-
# values corresponding to the domain between the two continuum maxima
82-
update_data[which((lambda_tmp - latest_intercept_tmp) >= 0 |
83-
latest_intercept_tmp == spectral_bands[nb_bands])] <- 0
84-
# values to eliminate for the next analysis: all spectral bands before latest intercept
85-
still_need_cr_tmp[which((lambda_tmp - latest_intercept_tmp) < 0)] <- 0
86-
# the max slope is known, as well as the coordinates of the beginning and ending
87-
# a matrix now has to be built
88-
convex_hull_tmp <- convex_hull_tmp +
89-
update_data * (intercept_hull_tmp + sweep((lambda_tmp - latest_intercept_tmp),
90-
1, slope[index_max_slope], "*"))
91-
# update variables
92-
convex_hull[sel,] <- convex_hull_tmp
93-
still_need_cr[sel,] <- still_need_cr_tmp
94-
lambda[sel,] <- lambda_tmp
95-
latest_intercept[sel,] <- latest_intercept_tmp
96-
intercept_hull[sel,] <- intercept_hull_tmp
97-
}
98-
cr_results0 <- mat_init[, 2:(nb_bands - 2)] / convex_hull[, 2:(nb_bands-2)]
99-
cr_results <- matrix(0, ncol = (nb_bands - 3), nrow = nb_samples)
100-
cr_results[cr_data$samples_to_keep, ] <- as.matrix(cr_results0)
26+
cr_results <- matrix(0, ncol = (nb_bands - 2), nrow = nb_samples)
27+
cr_results0 <- continuumRemoval(X = cr_data$mat_init,
28+
wav = cr_data$spectral_bands)
29+
cr_results[cr_data$samples_to_keep, ] <- as.matrix(cr_results0[,2:(nb_bands-1)])
30+
# # initialization:
31+
# # spectral band corresponding to each element of the data matrix
32+
# lambda <- repmat(matrix(spectral_bands, nrow = 1), nb_samples_update, 1)
33+
# # prepare matrices used to check evolution of the CR process
34+
# # - elements still not processed through continuum removal: initialization to 1
35+
# still_need_cr <- matrix(1, nrow = nb_samples_update, ncol = nb_bands)
36+
# # - value of the convex hull: initially set to 0
37+
# convex_hull <- matrix(0, nrow = nb_samples_update, ncol = nb_bands)
38+
# # - reflectance value for latest interception with convex hull:
39+
# # initialization to value of the first reflectance measurement
40+
# intercept_hull <- repmat(matrix(mat_init[, 1], ncol = 1), 1, nb_bands)
41+
# # - spectral band of latest interception
42+
# latest_intercept <- repmat(X = matrix(spectral_bands[1], ncol = 1),
43+
# m = nb_samples_update, n = nb_bands)
44+
# # number of spectral bands found as intercept
45+
# nb_intercept <- 0
46+
# # continues until arbitrary stopping criterion:
47+
# # stops when reach last spectral band (all values before last = 0)
48+
# # while (max(still_need_cr[, seq_len(nb_bands - 2)]) == 1 & (nb_intercept <= (nb_bands / 2))) {
49+
# while (max(still_need_cr[, seq_len((nb_bands - 2))]) == 1) {
50+
# nb_intercept <- nb_intercept + 1
51+
# # identify samples still needing continuum removal
52+
# sel <- which(still_need_cr[,(nb_bands-2)]==1)
53+
# # update variables to process samples needing CR only
54+
# nb_samples_update_tmp <- length(sel)
55+
# lambda_tmp <- lambda[sel,]
56+
# mat_init_tmp <- mat_init[sel,]
57+
# latest_intercept_tmp <- latest_intercept[sel,]
58+
# still_need_cr_tmp <- still_need_cr[sel,]
59+
# convex_hull_tmp <- convex_hull[sel,]
60+
# intercept_hull_tmp <- intercept_hull[sel,]
61+
# # Mstep give the position of the values to be updated
62+
# update_data <- matrix(1, nrow = nb_samples_update_tmp, ncol = nb_bands)
63+
# update_data[, nb_bands] <- 0
64+
# # initial step: first column set to 0; following steps: all bands below
65+
# # max of the convex hull are set to 0
66+
# update_data[which((lambda_tmp - latest_intercept_tmp) < 0)] <- 0
67+
# # compute slope for each coordinate
68+
# slope <- as.matrix((mat_init_tmp - intercept_hull_tmp) / (lambda_tmp - latest_intercept_tmp) * still_need_cr_tmp)
69+
# # set current spectral band and previous bands to -9999
70+
# if (!length(which(still_need_cr_tmp == 0)) == 0) {
71+
# slope[which(still_need_cr_tmp == 0)] <- -9999
72+
# }
73+
# if (!length(which(is.na(slope))) == 0) {
74+
# slope[which(is.na(slope))] <- -9999
75+
# }
76+
# # get max index for each row and convert into linear index
77+
# index_max_slope <- RowToLinear(max.col(slope, ties.method = "last"),
78+
# nb_samples_update_tmp, nb_bands)
79+
# # !!!! OPTIM: replace repmat with column operation
80+
# # update coordinates of latest intercept
81+
# latest_intercept_tmp <- repmat(matrix(lambda_tmp[index_max_slope],
82+
# ncol = 1), 1, nb_bands)
83+
# # update latest intercept
84+
# intercept_hull_tmp <- repmat(matrix(as.matrix(mat_init_tmp)[index_max_slope],
85+
# ncol = 1), 1, nb_bands)
86+
# # values corresponding to the domain between the two continuum maxima
87+
# update_data[which((lambda_tmp - latest_intercept_tmp) >= 0 |
88+
# latest_intercept_tmp == spectral_bands[nb_bands])] <- 0
89+
# # values to eliminate for the next analysis: all spectral bands before latest intercept
90+
# still_need_cr_tmp[which((lambda_tmp - latest_intercept_tmp) < 0)] <- 0
91+
# # the max slope is known, as well as the coordinates of the beginning and ending
92+
# # a matrix now has to be built
93+
# convex_hull_tmp <- convex_hull_tmp +
94+
# update_data * (intercept_hull_tmp + sweep((lambda_tmp - latest_intercept_tmp),
95+
# 1, slope[index_max_slope], "*"))
96+
# # update variables
97+
# convex_hull[sel,] <- convex_hull_tmp
98+
# still_need_cr[sel,] <- still_need_cr_tmp
99+
# lambda[sel,] <- lambda_tmp
100+
# latest_intercept[sel,] <- latest_intercept_tmp
101+
# intercept_hull[sel,] <- intercept_hull_tmp
102+
# }
103+
# cr_results0 <- mat_init[, 2:(nb_bands - 2)] / convex_hull[, 2:(nb_bands-2)]
104+
# cr_results <- matrix(0, ncol = (nb_bands - 3), nrow = nb_samples)
105+
# cr_results[cr_data$samples_to_keep, ] <- as.matrix(cr_results0)
101106
} else {
102107
cr_results <- matrix(0, ncol = (nb_bands - 3), nrow = nb_samples)
103108
}

‎R/crfun.R‎

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
#' computes continuum removal for individual spectra
2+
#' @description adapted from propsectr function
3+
#' https://github.com/l-ramirez-lopez/prospectr/blob/main/R/continuumRemoval.R
4+
5+
#' @param x numeric. original spectrum
6+
#' @param wav numeric. central wavelength for the spectral bands
7+
#' @param interpol character. linear or spline
8+
#
9+
#' @return continuum removed spectrum
10+
#' @importFrom grDevices chull
11+
#' @importFrom stats approx splinefun
12+
#' @export
13+
14+
crfun <- function(x, wav, interpol) {
15+
# need to define close neighbors corresponding to lower and upper bands
16+
neighbor <- wav[1]/1e4
17+
# convex hull
18+
id <- sort(chull(c(wav[1] - neighbor, wav, wav[length(wav)] + neighbor), c(-max(x), x, -max(x))))
19+
id <- id[-c(1, length(id))] - 1
20+
cont <- switch(interpol,
21+
linear = {
22+
approx(x = wav[id], y = x[id], xout = wav, method = "linear")$y
23+
},
24+
spline = {
25+
splinefun(x = wav[id], y = x[id])(wav)
26+
}
27+
)
28+
return(cont)
29+
}

‎R/filter_prior_CR.R‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -21,9 +21,9 @@ filter_prior_cr <- function(mat_init, spectral_bands) {
2121
samples_to_keep <- which(!null_sd == 0 & !is.na(null_sd))
2222
mat_init <- mat_init[samples_to_keep, ]
2323
# add negative values to the last column and update spectral bands
24-
mat_init <- cbind(mat_init, -9999)
24+
# mat_init <- cbind(mat_init, -9999)
25+
# spectral_bands <- c(spectral_bands, spectral_bands[nb_bands-1] + 100)
2526
nb_bands <- ncol(mat_init)
26-
spectral_bands <- c(spectral_bands, spectral_bands[nb_bands-1] + 100)
2727
return(list("mat_init" = mat_init, "spectral_bands" = spectral_bands,
2828
"nb_samples" = nb_samples, "samples_to_keep" = samples_to_keep))
2929
}

0 commit comments

Comments
 (0)