[R-sig-Geo] Needing to speed up a process involving calc() and cover() raster functions
Mathieu Rajerison
mathieu.rajerison at gmail.com
Wed Dec 23 15:04:47 CET 2015
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(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(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>:
> 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
>> 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
>
> _______________________________________________
> 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