 <p>Three steps:</p> <p>First, find all the cells in each polygon, return a list of 2-column matrices with the cell number and the value:</p> <pre><code>require(plyr) # for llply, laply in a bit... cell_value = extract(dat, polys,cellnumbers=TRUE) head(cell_value[[1]]) cell value [1,] 31 108 [2,] 32 108 [3,] 33 110 [4,] 92 110 [5,] 93 110 [6,] 94 111 </code></pre> <p>Second, turn into a list of similar matrices but add the x and y coords:</p> <pre><code>cell_value_xy = llply(cell_value, function(x)cbind(x,xyFromCell(dat,x[,"cell"]))) head(cell_value_xy[[1]]) cell value x y [1,] 31 108 8.581164 14.71973 [2,] 32 108 8.669893 14.71973 [3,] 33 110 8.758623 14.71973 [4,] 92 110 8.581164 14.67428 [5,] 93 110 8.669893 14.67428 [6,] 94 111 8.758623 14.67428 </code></pre> <p>Third, compute the weighted mean coordinate. This neglects any edge effects and assumes all grid cells are the same size:</p> <pre><code>centr = laply(cell_value_xy, function(m){c(weighted.mean(m[,3],m[,2]), weighted.mean(m[,4],m[,2]))}) head(centr) 1 2 [1,] 8.816277 14.35309 [2,] 8.327463 14.02354 [3,] 8.993655 13.82518 [4,] 8.467312 13.71929 [5,] 9.011808 13.28719 [6,] 9.745000 13.47444 </code></pre> <p>Now <code>centr</code> is a 2-column matrix. In your example its very close to <code>coordinates(polys)</code> so I'd make a contrived example with some extreme weights to make sure its working as expected.</p>