In R, how to average spatial points data over spatial grid squares

Question

Managed to solve problem now

I have a set of around 50 thousand points that have coordinates and one value associated with them. I would like to be able to place points into a grid averaging the associated value of all points that fall into a grid square. So I want to end up with an object that identifies each grid square and gives the average inside the grid square.

I have the data in a spatial points data frame and a spatial grid object if that helps.

Improving answer: I have definitely done some searching, sorry about the initial state of the question I had only managed to frame the question inside my own head; hadn't had to communicate it to anyone else before...

Here is example data that hopefully illustrates the problem more clearly

##make some data
longi <- runif(100,0,10)
lati <- runif(100,0,10)
value <- runif(500,20,30)

##put in data frame then change to spatial data frame
df <- data.frame("lon"=longi,"lat"=lati,"val"=value)
coordinates(df) <- c("lon","lat")
proj4string(df) <- CRS("+proj=longlat")

##create a grid that bounds the data
grd <- GridTopology(cellcentre.offset=bbox(df)[,1],
cellsize=c(1,1),cells.dim=c(11,11))
sg <- SpatialGrid(grd)

Then I hope to get an object albeit a vector/data frame/list that gives me the average of value in each grid cell/square and some way of identifying which cell it is.

Solution

##convert the grid into a polygon##
polys <- as.SpatialPolygons.GridTopology(grd) 
proj4string(polys) <- CRS("+proj=longlat")

##can now use the function over to select the correct points and average them
results <- rep(0, length(polys))

for(i in 1:length(polys)) {
  results[i] = mean(df$val[which(!is.na(over(x=df,y=polys[i])))])
}

My question now is if this is the best way to do it or is there a more efficient way?

Phrenic

Have you done any searching? I'm pretty sure you will find this is a duplicated question. If you do search and come up empty, then you should psot the search strategies that failed and then post a data example (not 50K points, maybe 50 x 50). – Kiwanis

Sorry, but this is too vague as is. For instance how do you define a "cell"? At the very least, you need to post sample data and show us what the desired output would look like. If you want to post the full dataset (fine by me), then upload it somewhere and post a link. – Barghest

Answer 1

Your description is vague at best. Please try to ask more specific answers preferably, with code illustrating what you have already tried. Averaging a single value in your point data or a single raster cell makes absolutely no sense.

The best guess at an answer I can provide is to use raster extract() to assign the raster values to a sp point object and then use tapply() to aggregate the values to your grouping values in the points. You can use the coordinates of the points to identify cell location or alternately, the cellnumbers returned from extract (per below example).

require(raster)
require(sp)

# Create example data
r <- raster(ncol=500, nrow=500)
  r[] <- runif(ncell(r))
    pts <- sampleRandom(r, 100, sp=TRUE)  

# Add a grouping value to points 
pts@data <- data.frame(ID=rownames(pts@data), group=c( rep(1,25),rep(2,25),
                       rep(3,25),rep(4,25)) )          

# Extract raster values and add to @data slot dataframe. Note, the "cells" 
#   attribute indicates the cell index in the raster. 
pts@data <- data.frame(pts@data, extract(r, pts, cellnumbers=TRUE))
  head(pts@data)

# Use tapply to cal group means  
tapply(pts@data$layer, pts@data$group, FUN=mean)

Hegumen

Thanks, I think I was wrong though in using the word cell, I mean the grid square. I have edited my question to try to clarify. I would like to group the points by which grid square they are in and average. – Phrenic

Just use the code I provided but aggregate by the cellnumber, which represents the unique ID of the raster cell. tapply(pts@data$MyValue, pts@data$cells, FUN=mean) – Hegumen