Needing to speed up a process involving calc() and cover() raster functions

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]]

system.time(fun1())

  system.time(fun2())

  system.time(fun3())

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]]

--
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

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]]

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]]