0

I want to select small region in the first raster and calculate the spatial average for that region and do the same for the other 11 rasters. so will finally get 12 values.

I tried this:

sami<- list.files("C:\\New folder (3)", "*.envi", full.names = TRUE)
saf=stack(sami)
plot(saf, 1)   ## to select a region
e <- drawExtent()## I selected the region
saf_test <- crop(saf, e)

is this right to do so?

Then how I calculate the spatial average for the selected region?and do the same for all other rasters. Thanks in advance.

hyat
  • 1,047
  • 4
  • 15
  • 34
  • Please make your situation reproducible, i.e. provide us with the data and the code needed to mimic your situation. See http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example for more tips on how to do this. – Paul Hiemstra Mar 20 '13 at 11:11
  • did you try Zonal statistics in Raster Library? zonal(x, z, stat=’mean’, digits=0, na.rm=TRUE, ...) – Gianni Spear Mar 20 '13 at 11:14
  • What about `stackApply` You could use `extent` to select the region. This will make it a bit more reproducible (and easier in the long run). You can check the raster vignette for some tips on how to generate some sample data to make your question reproducible. – Roman Luštrik Mar 20 '13 at 11:16
  • @ZadSim it's not important your real data but just a code as r <- raster(ncols=10, nrows=10), etc etc – Gianni Spear Mar 20 '13 at 11:31
  • @ZadSim In order to stop attracting downvotes to your questions (I don't think you have been asking *bad* questions) you need to focus on making your questions reproducible as Paul points out above. Simply adding a link to you data is not good enough - it forces people to download it manually. Better to create some code & dummy data that recreates your problem in R so people can copy paste the code from your question and develop it. Check the docs of the package you are using, there is usually some sample data or methods to create sample data (especially in the excellent raster package). – Simon O'Hanlon Mar 20 '13 at 11:51

1 Answers1

1

You can do it like this. I use some example data from the raster package. Try it and adapt it to your raster stack (this should be easy). The extract function will return a matrix. One column of values for the region selected by drawExtent for each layer, and you can then use colMeans for the mean of that area:

    #This particular stack has 3 'layers' - one each for the red, green and blue channels of the picture
    saf <- stack(system.file("external/rlogo.grd", package="raster")) 
    plotRGB( saf )
    e <- drawExtent()

enter image description here

    vals <- extract( saf , e , nl = nlayers( saf ) )
    vals <- colMeans( vals )
    vals
#   red    green     blue 
#   185.9368 191.9158 208.7825
Simon O'Hanlon
  • 58,647
  • 14
  • 142
  • 184
  • You can manually create an extent object. If you draw the extent, then type `e` into the console (or whatever you called your extent) it will bring up the values of your extent. Now, in your code use those values to create an extent like so... `e <- extent(xmin, xmax, ymin, ymax)` inputting the correct values in the specified order. Does that make sense? – Simon O'Hanlon Mar 20 '13 at 12:33
  • @ZadSim good grief, always more questions and still no accept or upvote?! `plot( saf )` , `plot( e , add = TRUE ) ` – Simon O'Hanlon Mar 20 '13 at 12:42
  • @ZadSim thanks! :-), also to get the exact value of the extent you may want to do `e <- drawExtent()` , then `paste( e@xmin )` , `paste( e@xmax )` etc. – Simon O'Hanlon Mar 20 '13 at 12:45
  • @ZadSim what was the error? That should have worked. Perhaps you have some columns which are all NA? – Simon O'Hanlon Mar 20 '13 at 12:57
  • Are you sure you are running colMeans on a matrix, and not on the colMeans of a previous calculation?!! if you type in `vals` do you get a mtrix or a vector of values? Try changing the code to `vals <- extract(...)` , `valMean <- colMeans(...)` to avoid any confusion and run it again from the beginning! – Simon O'Hanlon Mar 20 '13 at 13:16