I've been working with some raster data (19 years of precipitation data) and a shapefile. The raster extent represents a square-shaped area and the shapefile is within the extent of the raster.

enter image description here

I want to obtain the average monthly rainfall, I mean the average rainfall for each month for all the time series.

Each month has 19 relevant rasters (from the time series period) and I can get the average rainfall which results to be a single raster for each month but don't know how to get the areal average rainfall (a single value that represents the rainfall within the shapefile).

I loaded the rasters as single objects in R. The raster function mean(r1, r2, r3) gives a raster as output, I want a single value that represents only the area of study.

So far I've been working with R and was wondering if I can get that calculation using the raster package or any other package.

  • 1
    Define your extent equal to one of rectangle and run cell statistics on relevant rasters
    – FelixIP
    Jul 15, 2017 at 0:57
  • Raster statistics on output
    – FelixIP
    Jul 15, 2017 at 2:58
  • Have you read the rasters in? Are they a number of single-band raster objects or have you read them in as a multi-band stack object? Did you know you can just do mean(r1,r2,r3) to get the raster of cell-means of a number of rasters?
    – Spacedman
    Jul 15, 2017 at 8:05
  • 1
    I know, thank you for answering. I added more details. The mean() function works for the entire extent of the raster and gives a raster as output. I need a single value representing the inside area of the shapefile.
    – noriega
    Jul 15, 2017 at 13:53
  • extract(stack(raster1, raster2, ...), polygon, fun= mean, na.rm= TRUE) perhaps?
    – Kazuhito
    Jul 16, 2017 at 8:24

1 Answer 1


Thanks to the comments I managed to answer my own question.

I've been working these past days and the code shared by Kazuhito in the comments:

extract(stack(raster1, raster2, ...), polygon, fun= mean, na.rm= TRUE)

happens to work perfectly and thereby is the answer I was looking for.

I just want to point out some details. Despite the code is very simple and straightforward, in order to convert a shapefile into a polygon I did this:

kis_stk = stack(raster1, raster2, ....)
rgdal::readOGR('directory-of-shapefile', stringsAsFactors = F) %>%
                sp::spTransform(CRS(projection(kis_stk[[1]]))) -> basin
extract(kis_stk, spPolygons(basin), fun= mean, na.rm= TRUE)

This worked with a .shp file I had in the directory-of-shapefile location.

Hope all this info helps someone.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.