[Gdal-dev] mean value of rasters

[Roger Bivand]

On Fri, 5 May 2006, Christopher Condit wrote:

> I tried in R and came very close.  My script looked something like this
> (just to make sure things were working):
> 
> library(rgdal)
> ds <- GDAL.open("C:/path/test.tif")
> vals <- getRasterData(ds)
> GDAL.close(ds)
> write.table(vals, "C:/path/out.txt")
> 
> But when I look at that file as an ASCII grid (which makes sense to do,
> right?) my original image is rotated 90 degrees and doubled.
> Am I missing something about the rgdal component?  Seems that it should
> have worked...

Sense depends on which way you look at things ...

Look at for example:

library(rgdal)
cea <- system.file("pictures/cea.tif", package = "rgdal")[1]
system(paste("gdalinfo -mm", cea))
cea_SGDF <- readGDAL(cea)
summary(cea_SGDF)
image(cea_SGDF)

where cea_SGDF is a SpatialGridDataFrame with slots:

> getSlots(class(cea_SGDF))
           data            grid      grid.index          coords          bbox 
"AttributeList"  "GridTopology"       "integer"        "matrix"      "matrix" 
    proj4string 
          "CRS" 

If you have N images of an identical region, the images are so large that
you don't want all in memory at once, and you have full control of the
NAs, you can do something like:

cumulator <- readGDAL("image1.tif")
for (i in 2:N) {
  tmpimage <- readGDAL(paste("image", i, ".tif", sep=""))
  cumulator$band1 <- cumulator$band1 + tmpimage$band1
}
cumulator$band1 <- cumulator$band1/N
writeGDAL(cumulator, "cumulator.tif")

But this is fragile - any NAs wil NA the whole cell, and if the lengths of
band1 differ from image to image, all bets are off (you say that the
region is the same for all images, hence the risky suggestion above). The
GRASS route provides GIS-style security if you need it, if you don't,
check the data within R first, possibly on a sub-region (offset and
region.dim arguments to readGDAL), to find out what is going on - maybe a 
print(summary(tmpimage)) in the loop.

Roger

> 
> Thanks,
> 
> Chris
> 
> 
> -----Original Message-----
> From: Frank Warmerdam [mailto:fwarmerdam at gmail.com] On Behalf Of Frank
> Warmerdam
> Sent: Thursday, May 04, 2006 1:00 PM
> To: Christopher Condit
> Cc: Gdal-dev at lists.maptools.org
> Subject: Re: [Gdal-dev] mean value of rasters
> 
> Christopher Condit wrote:
> > I've got N raster files (GeoTiffs) of equal size and dimensions.  I'd
> > like to create a single raster with the mean value of these N rasters
> > (cell by cell).  Is there a good way to do this with GDAL, or with
> > another tool?  
> 
> Chris,
> 
> There is not direct way of doing this with the commandline utilities.
> It would be relatively easy to implement using Python bindings for
> GDAL and Numerical Python.  I imagine this could also be done in
> something
> like R or GRASS relatively directly.
> 
> Best regards,
> 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no