[R-sig-Geo] Needing to speed up a process involving calc() and cover() raster functions

Mathieu Rajerison mathieu.rajerison at gmail.com
Mon Jan 4 11:51:13 CET 2016


Hi,

And thanks for all the answers.

I've tested Robert's elegant solution successfully.
I haven't yet tested your solution, Banjamin, but it seems a really good
improvement for speeding all this up and I'll soon give it a try. Thanks
for your awesome RSToolbox package

To answer your question "why setting all pixels to 1", it's just that
afterwards, I smooth the result with a mean-shift focal function. This way,
it gives all the extracted pixels the same value/weight.

Thanks again and happy new yeaR

Mathieu


2015-12-25 6:46 GMT+01:00 Robert J. Hijmans <r.hijmans at gmail.com>:

> How about this:
>
> b <- brick(system.file("external/rlogo.grd", package="raster"))
> classified <- b / 255
>
> threshs <- c(.1, .3, .5)
> x <- classified < threshs
> y <- max(x, na.rm=TRUE)
> z <- reclassify(y, cbind(0, NA))
>
>
> Robert
>
>
> On Tue, Dec 22, 2015 at 1:57 AM, Mathieu Rajerison
> <mathieu.rajerison at gmail.com> 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
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

	[[alternative HTML version deleted]]



More information about the R-sig-Geo mailing list