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