# Calculating average monthly rainfall from raster using shapefile

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.

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.

• Define your extent equal to one of rectangle and run cell statistics on relevant rasters Jul 15, 2017 at 0:57
• Raster statistics on output 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? Jul 15, 2017 at 8:05
• 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. Jul 15, 2017 at 13:53
• `extract(stack(raster1, raster2, ...), polygon, fun= mean, na.rm= TRUE)` perhaps? Jul 16, 2017 at 8:24

## 1 Answer

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.