sum should be faster than calc(x, sum) because sum uses rowSums whereas calc, in this case, would call apply(values(x), 1, sum). calc(x, function(i) rowSums(i)) should be equivalent to sum (in speed and results).
Whether making a RasterBrick is worth it depends on how many times you will use these files for this type of computation. If it is only for one computation it should be slower to first make bricks, as the main bottleneck in these computations is disk I/O. I cannot think of another way to speed things up right now --- something we are working through clustering (probably not too helpful here), and C routines. 30 hours is long, but it is the computer that's working, not you. Robert On Sun, May 13, 2012 at 4:46 AM, Matteo Mattiuzzi < matteo.mattiu...@boku.ac.at> wrote: > Dear Thiago, > > I don't know what makes sum() more suitable comparing to calc() but I > think a fakt that you are useing a stack (and not a brick) influences also > the speed of processing: > > ###################### > library(raster) > r <- raster() > r[] <- 1:ncell(r) > > setOptions(todisk=TRUE) > system.time(st <- stack(r,r,r,r,r,r,r,r,r)) > system.time(br <- brick(r,r,r,r,r,r,r,r,r)) > > system.time(sum(st,na.rm=FALSE)) > system.time(sum(br,na.rm=FALSE)) > > system.time(calc(st,sum,na.rm=FALSE)) > system.time(calc(br,sum,na.rm=FALSE)) > > setOptions(todisk=FALSE) > ####################### > > # Of coarse the creation of a brick takes a bit of time since the data > needs to be re-organised physically on the disk. > # But this initial overhead worthwhile relatively quickly, you have to > check this. > > # Cheers Matteo > > > > > >>> Thiago Veloso 13.05.12 1.53 Uhr >>> > Robert, > Thanks for your input. Your suggestion is more suitable to my case, but > it takes almost one hour to sum a stack containing 46 satellite images. > Please see below: > > sclass : RasterStack dimensions : 5568, 8289, 46153152, 46 > (nrow, ncol, ncell, nlayers)resolution : 0.00898, 0.00898 (x, y)extent > : -104.4326, -29.99736, -40.00064, 10 (xmin, xmax, ymin, ymax)coord. > ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 min > values : -32768, -32768, -32768, -32768, -32768, -32768, -32768, -32768, > -32768, -32768, -32768, -32768, -32768, -32768, -32768, ... max values : > 32767, 32767, 32767, 32767, 32767, 32767, 32767, 32767, 32767, 32767, > 32767, 32767, 32767, 32767, 32767, ... layer names : gpp2004001.Gpp_1km, > gpp2004009.Gpp_1km, gpp2004017.Gpp_1km, gpp2004025.Gpp_1km, > gpp2004033.Gpp_1km, gpp2004041.Gpp_1km, gpp2004049.Gpp_1km, > gpp2004057.Gpp_1km, gpp2004065.Gpp_1km, gpp2004073.Gpp_1km, > gpp2004081.Gpp_1km, gpp2004089.Gpp_1km, gpp2004097.Gpp_1km, > gpp2004105.Gpp_1km, gpp2004113.Gpp_1km, ... > For reference, total size of the 46 images which compose the stack > reaches 4GB. Here is how long it takes to perform the sum over the stack: > > system.time(sos <- sum(s,na.rm=FALSE)) user system elapsed > 743.098 31.258 3104.113 > The machine I am using to process this dataset is an 3GHz quadcore with > 16GB of RAM, running Opensuse and latest R (2.15.0). > As processing my entire data involves performing this operation at least > 30 times, is there a way to optimize the process? > All the best, Thiago. > --- On Sat, 12/5/12, Robert J. Hijmans wrote: > > From: Robert J. Hijmans > Subject: Re: [R-sig-Geo] Antw: Error when using rowSums along with a > RasterStack object > To: thi_vel...@yahoo.com.br > Cc: r-sig-geo@r-project.org > Date: Saturday, 12 May, 2012, 13:43 > > You can also do > sos <- sum(s, na.rm=FALSE) > Best, Robert > > > On Sat, May 12, 2012 at 2:26 AM, Matteo Mattiuzzi wrote: > > Dear Thiago, > > > > is this what you need? > > #### > > library(raster) > > r <- raster() > > r[]<-1:ncell(r) > > r2 <- raster() > > r2[]<-ncell(r):1 > > r <-brick(r,r2) > > plot(r) > > > > sos <- calc(r,sum) > > x11() > > plot(sos) > > #### > > Matteo > > > > sos <- calc(r,sum) > > plot(sos) > > > > >>> Thiago Veloso 11.05.12 23.29 Uhr >>> > > Dear R-colleagues, > > > > I am trying to sum all layers of a RasterStack object, which is summarized > below: > > > > > s > > class : RasterStack > > dimensions : 5568, 8289, 46153152, 46 (nrow, ncol, ncell, nlayers) > > resolution : 0.00898, 0.00898 (x, y) > > extent : -104.4326, -29.99736, -40.00064, 10 (xmin, xmax, ymin, ymax) > > coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 > +towgs84=0,0,0 > > min values : -32768, -32768, -32768, -32768, -32768, -32768, -32768, > -32768, -32768, -32768, -32768, -32768, -32768, -32768, -32768, ... > > max values : 32767, 32767, 32767, 32767, 32767, 32767, 32767, 32767, > 32767, 32767, 32767, 32767, 32767, 32767, 32767, ... > > layer names : gpp2004001.Gpp_1km, gpp2004009.Gpp_1km, gpp2004017.Gpp_1km, > gpp2004025.Gpp_1km, gpp2004033.Gpp_1km, gpp2004041.Gpp_1km, > gpp2004049.Gpp_1km, gpp2004057.Gpp_1km, gpp2004065.Gpp_1km, > gpp2004073.Gpp_1km, gpp2004081.Gpp_1km, gpp2004089.Gpp_1km, > gpp2004097.Gpp_1km, gpp2004105.Gpp_1km, gpp2004113.Gpp_1km, ... > > > > > The raster 's' is composed by 46 satellite images at a resolution of 1km > (0.00898degree). Each layer is a 8-day composition of leaf area index and I > need to integrate an entire year of LAI. When I try to run the "rowSums" > function in this RasterStack, the following error is displayed: > > > > > > sos <- rowSums(s,na.rm=FALSE) > > Error in rowSums(s, na.rm = FALSE) : > > 'x' must be an array of at least two dimensions > > > > The same occurs when I try to use .rowSums providing matrix dimensions: > > > > > sos <- .rowSums(s,5568,8289) > > Error in .rowSums(s, 5568, 8289) : 'x' must be numeric > > > > > > What am I doing wrong? Is there another way to integrate the layers of a > RasterStack object? > > > > Thanks in advance, > > Thiago Veloso. > > > > P.S.: Attaching sessionInfo output to provide useful information: > > > > > sessionInfo() > > R version 2.15.0 (2012-03-30) > > Platform: x86_64-unknown-linux-gnu (64-bit) > > > > locale: > > [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C > > [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 > > [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 > > [7] LC_PAPER=C LC_NAME=C > > [9] LC_ADDRESS=C LC_TELEPHONE=C > > [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C > > > > attached base packages: > > [1] stats graphics grDevices utils datasets methods base > > > > other attached packages: > > [1] ncdf_1.6.6 rgdal_0.7-11 raster_1.9-92 sp_0.9-99 > > > > loaded via a namespace (and not attached): > > [1] grid_2.15.0 lattice_0.20-6 > > > > > > _______________________________________________ > > R-sig-Geo mailing list > > R-sig-Geo@r-project.org > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > > > > > > [[alternative HTML version deleted]] > > > > _______________________________________________ > > R-sig-Geo mailing list > > R-sig-Geo@r-project.org > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > > > > > [[alternative HTML version deleted]] > > > > [[alternative HTML version deleted]] > > _______________________________________________ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
