[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