[R-sig-Geo] Rasters average

[bart]

I think the calculations could be a lot quicker if you make a raster 
stack or raster brick from your seven rasters and use that. see example 
below. Maybe you have to play around with your in memory settings and 
use a file on the disk. But i'm still able to do it using rasters upto 
200000 points easily in memory within seconds

 > (a<-brick(raster(matrix(runif(10), ncol=2)), 
raster(matrix(c(runif(8), NA, NA), ncol=2))))
class       : RasterBrick
dimensions  : 5, 2, 2  (nrow, ncol, nlayers)
ncell       : 10
resolution  : 0.5, 0.2  (x, y)
projection  : NA
extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
values      : from memory
min values  : 0.015 0.015
max values  : 0.94 0.95

 > values(a[[2]])
  [1] 0.13370844 0.94727248 0.45225123 0.16531864 0.62212769 0.01502728
  [7] 0.66552310         NA 0.72266646         NA
 > values(a[[1]])
  [1] 0.41511274 0.90717270 0.59899444 0.17772793 0.93664444 0.53157684
  [7] 0.05073588 0.91242408 0.01484997 0.29561654
 > values(mean(a, na.rm=TRUE))
  [1] 0.2744106 0.9272226 0.5256228 0.1715233 0.7793861 0.2733021 0.3581295
  [8] 0.9124241 0.3687582 0.2956165
 >


On 02/09/2011 11:05 AM, luca candeloro wrote:
> Hi all,
> I'm trying to use MODIS 1day raster (.tif file) to obtain a, say, weekly
> raster which is the average of the seven 1 day raster.
> Using the following code (library raster) I can't obtain the results due to
> computational limits.
>
> #Input files (same dimension: [1310,1391]
>
> Map.raster1=raster("C:/MOD11A2/Confronto/1d001.tif")
> Map.raster2=raster("C:/MOD11A2/Confronto/1d002.tif")
> Map.raster3=raster("C:/MOD11A2/Confronto/1d003.tif")
> Map.raster4=raster("C:/MOD11A2/Confronto/1d004.tif")
> Map.raster5=raster("C:/MOD11A2/Confronto/1d005.tif")
> Map.raster6=raster("C:/MOD11A2/Confronto/1d006.tif")
> Map.raster7=raster("C:/MOD11A2/Confronto/1d007.tif")
>
>
> for (r in 1:1310) {
> for (c in 1:1391) {
> if(Map.raster1[r,c]==0) v1=NA else v1= Map.raster1[r,c]
> if(Map.raster2[r,c]==0) v2=NA else v2= Map.raster2[r,c]
> if(Map.raster3[r,c]==0) v3=NA else v3= Map.raster3[r,c]
> if(Map.raster4[r,c]==0) v4=NA else v4= Map.raster4[r,c]
> if(Map.raster5[r,c]==0) v5=NA else v5= Map.raster5[r,c]
> if(Map.raster6[r,c]==0) v6=NA else v6= Map.raster6[r,c]
> if(Map.raster7[r,c]==0) v7=NA else v7= Map.raster7[r,c]
> media=mean(c(v1,v2,v3,v4,v5,v6,v7),na.rm=TRUE)
> if (is.na(media)) media=0
> Map.raster[r,c]=media
> }
> }
>
> Is there an easiest way to average the values of "k" raster?
> Thanks,
> Luca.
>
> 	[[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
>