[R-sig-Geo] Needing to speed up a process involving calc() and cover() raster functions
Benjamin Leutner
benjamin.leutner at uni-wuerzburg.de
Wed Dec 23 18:04:00 CET 2015
so, essentially what you want is pixels which have NAs throughout all
layers (i.e. are above the thresholds) set to NA and the rest set to 1.
although, I don't quite see the use-case, but this is how you could do
it a little faster:
library(Rcpp) ## You need RcppArmadillo also installed
cppFunction("
arma::vec allNA(arma::mat x, arma::vec tr){
arma::vec out(x.n_rows);
out.fill(NA_REAL);
for(int i = 0; i < tr.n_elem; i++) {
out.elem(find(x.col(i) < tr(i))).ones();
}
return out;
}
", depends = "RcppArmadillo")
r <- calc(classified, function(x){allNA(x,tr=threshs)}, forcefun=T)
-
Benjamin
On 23.12.2015 15:04, Mathieu Rajerison wrote:
> Hi Benjamin,
>
> Thanks for your answer
>
> Sorry for not being precise. Actually, I want to threshold each band
> given thresholds. I give the 1 value to each pixel, in each band, that
> is below the threshold value.
>
> Then, I want to cover all these bands into one, in order to combine
> all the layers.
>
> I tried different functions. Here are the different performances, if
> you want to have a look. Te second function below is the fastest one.
> If you think you have better, do not hesitate.
>
> ## FIRST FUNCTION
>
> fun1 = function() {
> threshs = c(.1,.2,.1)
> tamis = raster(classified); tamis[] = 0
> out = stack(lapply(1:nlayers(classified), function(i)
> cover(clamp(classified[[i]], upper = threshs[[i]], useValues = FALSE),
> tamis)))
> r = calc(out, sum)
> r[r[]==0] = NA; r[which(!is.na <http://is.na>(r[]))] = 1
> }
>
> ## SECOND
>
> fun2 = function() {
> threshs = c(.1,.2,.1)
> out = lapply(1:nlayers(classified), function(i)
> clamp(classified[[i]], upper = threshs[[i]], useValues = FALSE))
> r = out[[1]]
> for(i in 2:length(out)) {
> r = cover(out[[i]], r)
> }
> r[which(!is.na <http://is.na>(r[]))] = 1
> }
> ## THIRD
>
> fun3 = function() {
> threshs = c(.1,.2,.1)
> out=lapply(1:nlayers(classified), function(i) {
> out[[i]] = calc(classified[[i]], function(x) {x[x >= threshs[i]]
> = NA;
> x[x < threshs[i]] = 1;
> return(x)})
> })
> r = out[[1]]
> for(i in 2:length(out)) {
> r = cover(out[[i]], r)
> }
> }
> system.time(fun1())
> system.time(fun2())
> system.time(fun3())
>
> > system.time(fun1())
> user system elapsed
> 63.05 2.67 65.81
> > system.time(fun2())
> user system elapsed
> 43.14 0.88 *44.52 *
> > system.time(fun3())
> user system elapsed
> 48.67 1.51 50.62
>
>
> Best,
>
> Mathieu
>
>
>
> 2015-12-22 14:09 GMT+01:00 Benjamin Leutner
> <benjamin.leutner at uni-wuerzburg.de
> <mailto:benjamin.leutner at uni-wuerzburg.de>>:
>
> Hi Mathieu,
>
> your question is rather difficult to understand. From the context
> I gather that you are referring to the results of the sam()
> function from RStoolbox.
> Further, I assume you want to threshold each layer for a maximum
> spectral angle and then find the class with the minimum spectral
> angle per pixel, right?
>
> In this case you could do:
>
> out <- stack(lapply(1:nlayers(classified), function(i)
> clamp(classified[[i]], upper = threshs[[i]], useValues = FALSE)))
> class <- which.min(out)
>
> Cheers,
> Benjamin
>
>
> On 22.12.2015 10:57, Mathieu Rajerison wrote:
>
> Hi,
>
>
> I use RSToolBox to classify a RGB raster.
>
> I have a resulting RasterBrick which has as many layer as end
> members, in
> my case 3 for different tones of blue.
>
> I reclassify each band with calc to extract the pixels which
> have a small
> angle mapping value. The threshold used is different depending
> on the
> endmember layer.
>
> I finally assembly all the bands with the cover function.
>
> I needed to increase the memory limit assigned to R to have it
> worked. I
> suspect that my code could be optimized, but I don't know in
> which way.
>
> Here is the part of my code, that I think, could be
> optilmized, if you want
> to have a look and give some advice :
>
> # RECLASSIFY
> # classified is the classified RasterStack
> # here I change the values of each band to 1 or NA depending
> on the
> spectral angle mapping value.
> # Is calc() slower than reclassify() for this purpose as I
> have only one
> threshold value ?
>
> threshs = c(.1,.2,.1)
> for (i in 1:nlayers(classified )) {
>
> clas = classified[[i]]
> thresh=threshs[i]
>
> out[[i]] = calc(clas, function(x) {x[x >= thresh] = NA;
> x[x < thresh] = 1;
> return(x)})
> }
>
> # COVERING
> r = out[[1]]
> for(i in 2:length(out)) {
> r = cover(out[[i]], r) ## I cover by iteration
> }
>
> plot(r) # r is the final combined raster
>
> ================================================================
> My sessionInfo() :
>
> sessionInfo()
>
> R version 3.1.2 (2014-10-31)
> Platform: x86_64-w64-mingw32/x64 (64-bit)
>
> locale:
> [1] LC_COLLATE=French_France.1252 LC_CTYPE=French_France.1252
> LC_MONETARY=French_France.1252
> [4] LC_NUMERIC=C LC_TIME=French_France.1252
>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods
> base
>
> other attached packages:
> [1] R.utils_2.1.0 R.oo_1.19.0 R.methodsS3_1.7.0
> igraph_1.0.1 scatterplot3d_0.3-36
> [6] gdalUtils_2.0.1.7 spdep_0.5-88 Matrix_1.2-2
> maptools_0.8-34 spgrass6_0.8-8
> [11] XML_3.98-1.3 rgeos_0.3-8 FNN_1.1
> rgdal_0.9-2 RStoolbox_0.1.1
> [16] raster_2.3-40 sp_1.0-17
>
> loaded via a namespace (and not attached):
> [1] boot_1.3-13 car_2.0-25 caret_6.0-57
> coda_0.18-1
> codetools_0.2-9 colorspace_1.2-6
> [7] deldir_0.1-9 digest_0.6.8 doParallel_1.0.8
> foreach_1.4.2
> foreign_0.8-61 geosphere_1.4-3
> [13] ggplot2_1.0.1 grid_3.1.2 gtable_0.1.2
> iterators_1.0.7 lattice_0.20-29 LearnBayes_2.15
> [19] lme4_1.1-10 magrittr_1.5 MASS_7.3-35
> MatrixModels_0.4-1 mgcv_1.8-3 minqa_1.2.4
> [25] munsell_0.4.2 nlme_3.1-118 nloptr_1.0.4
> nnet_7.3-8
> parallel_3.1.2 pbkrtest_0.4-2
> [31] plyr_1.8.3 proto_0.3-10 quantreg_5.19
> Rcpp_0.12.0
> reshape2_1.4.1 scales_0.2.5
> [37] SparseM_1.7 splines_3.1.2 stats4_3.1.2
> stringi_0.5-5
> stringr_1.0.0 tools_3.1.2
>
> ============================================================
> Best,
>
> Mathieu
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org <mailto:R-sig-Geo at r-project.org>
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
> --
> Benjamin Leutner M.Sc.
>
> Department of Remote Sensing
> University of Wuerzburg
> Campus Hubland Nord 86
> 97074 Wuerzburg, Germany
>
> Tel: +49-(0)931-31 89594 <tel:%2B49-%280%29931-31%2089594>
> Fax: +49-(0)931-31 89594-0 <tel:%2B49-%280%29931-31%2089594-0>
> Email: benjamin.leutner at uni-wuerzburg.de
> <mailto:benjamin.leutner at uni-wuerzburg.de>
> Web: http://www.fernerkundung.uni-wuerzburg.de
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org <mailto:R-sig-Geo at r-project.org>
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
--
Benjamin Leutner M.Sc.
Department of Remote Sensing
University of Wuerzburg
Campus Hubland Nord 86
97074 Wuerzburg, Germany
Tel: +49-(0)931-31 89594
Fax: +49-(0)931-31 89594-0
Email: benjamin.leutner at uni-wuerzburg.de
Web: http://www.fernerkundung.uni-wuerzburg.de
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list