# Weighted mean of raster by polygons in R?

I have two rasters, raster_a and raster_b, and a shapefile, world_shapefile. Both rasters are of the world extent and overlay each other.

What I want to do is to create a weighted average of raster_a, where the weights are from raster_b, by each polygon in my world_shapefile.

I guess it is quite straightforward and goes something like: read in shapefile and rasters, create the weighed mean by polygons.

EDIT: My code currently looks as follows:

``````require(raster)
require(sp)
require(gdal)

setwd("C:/Users/XXX/Desktop/data"

polys <- readOGR(dsn=get_wd(), layer = "world_shapefile"
raster.a <- raster(raster_a.txt)
raster.b <- raster(raster_b.txt)

plot(polys)
plot(raster.a)
plot(raster.b)
``````

Everything looks perfect up to this point and then I get stuck...

• What does your code so far look like? – PolyGeo Jun 2 at 12:29
• Have you managed to read in the shapefile and the rasters? Show us the code, so at least we know where to start from. Better still create some sample data we can all use to demonstrate the solution. – Spacedman Jun 2 at 12:47
• I have updated my post - thanks! – Picko90 Jun 3 at 9:22
• Are your expected results a new raster or the weighted mean values for each polygon? – Jeffrey Evans Jun 3 at 15:26
• The output should be the weighted mean values for each polygon. – Picko90 Jun 3 at 15:48

You are barely at the start of this analysis. As such, I would highly recommend trying to track down any resources as a starting point to learn the basics of spatial analysis in R. Honestly, Google is a good start and it is clear that you did not search this StackExchange very well because the fundamentals of your task have been covered numerous times herein.

I will give you a general workflow, with some potential functions, and then you can ask for additional assistance once you get to a point where you have worked through some actual code.

1. Read rasters using the `raster` function in the raster package and polygon vector data with the `readOGR` function from rgdal (among other approaches that results in an SpatialPolygonsDataFrame object), validate projections.

2. Using the `extract` function in the raster package, extract raster values associated with each polygon.

3. Now, here it gets a bit more complicated. Multiple raster cells will intersect each polygon. As such, the results of the extract function are stored in a list object (as vectors or, in the case of raster stacks, data.frames). This data needs to be aggregated in some way so that a single value is represented. This is where your weighting comes in. Since you need a weighted mean, based on a second raster, you will need to get tricky with either a `for` loop with an index or even better, the `mapply` function that allows you to use an `lapply` (list apply) type function across multiple list objects.

I should also add that, depending on how you create your raster object, this could be simplified. If your rasters share common origin coordinates, resolution and rows/columns then you could read the data as a stack (more than one raster in a single object). In this case the results of the extract function will be a data.frame, with columns for each raster, for each polygon in a list. Accordingly, you can then just use lapply and pass a function that calls the correct column for the values and weights.

Something along the lines of:

``````r <- stack(raster_a.txt, raster_b.txt)
rv <- extract(r, polys)
unlist(lapply(rv, FUN = function(x) { weighted.mean(x=x[,1], w=x[,2], na.rm = TRUE) }))
``````

Unless your weight raster already represents weights, you will likely need to use an `lapply` function to create a new list object where weights to each polygon sum to one, the assumption of the `weighted.mean` function.

• Jeffrey, thanks for this. The idea is to somehow apply weights to each grid cell, in raster_a, by using the cosine of its latitude, which I believe is contained in raster_b. All I want to do is have have the average value of raster_a in country X whilst account for this cell area. Cheers. – Picko90 Jun 3 at 16:05