An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20110209/4bc24f76/attachment.pl>
Rasters average
4 messages · luca candeloro, Kranstauber, Bart, Tomislav Hengl +1 more
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
Try also using RSAGA. I usually process 100 millions of pixels without a problem. See: > rsaga.get.usage(lib="geostatistics_grid", module=5) module name: Statistics for Grids Here is an example: [http://spatial-analyst.net/scripts/getMOD12C1.R] HTH, T. Hengl http://www.wewur.wur.nl/popups/vcard.aspx?id=HENGL001 Op 9-2-2011 11:05, luca candeloro schreef:
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
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.
Luca,
You can do something like:
s <- stack(paste("C:/MOD11A2/Confronto/1d00", 1:7, ".tif", sep=""))
fun <- function(x, ...) { x[x==0] <- NA; mean(x, na.rm=TRUE) }
m <- calc(s, fun)
I am not sure what you mean by "computational limits"; are you saying that
this takes a long time? If so, that would be because you are looping over
cells. Avoid that, and things will be faster.
Robert
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Rasters-average-tp6007034p6008219.html Sent from the R-sig-geo mailing list archive at Nabble.com.