Hi Gavin,

I have been working with some code that does something very similar.  In my
case, I am intersecting points with grids of varying resolutions in order
to collect the grid identifier (and coordinates for grid cell centers),
then averaging the results of several variables held by the
SpatialPointsDataFrame by year (the points contain observations from
different years), then by grid identifier.  This relies on
doBy::summaryBy().

I am sorry that I don't have a a dummy set of points for you to run it
with, but perhaps there is something in here that helps you with what you
are trying to do.

Cheers, Lyndon

# Create sampling grid
# Grid template, bumped out 20 km on either side
rex <- extent(p)
gr <- raster(xmn = rex@xmin - 20000, xmx = rex@xmax + 20000, ymn =
rex@ymin- 20000, ymx = rex@ymax+ 20000)

# Adjust grid size and collecting grid ids under points for each resolution
polID.l <- c("km1", "km5", "km10", "km20")
res.l <- c(1000, 5000, 10000, 20000)
coord.l <- c("xy1", "xy5", "xy10", "xy20")
grid.l <- vector("list", length(res.l))  # list to hold dummy grids

# Extract grid cell IDs
for(i in 1:4) {
  gr.2 <- gr
  res(gr.2) <- res.l[i]
  print(res.l[i])
  projection(gr.2) <- p@proj4string
  gr.2[] <- 1:ncell(gr.2)
  grid.l[[i]] <- gr.2
  xy <- cbind(1:ncell(gr.2), xyFromCell(gr.2, 1:ncell(gr.2)))  # Center
coordinates
  colnames(xy)[1] <- polID.l[i]
  print(coord.l[i])
  assign(coord.l[i], xy)
  p$km <- extract(gr.2, p)
  colnames(p@data)[length(colnames(p@data))] <- polID.l[i]
}

# Create lists for lapply: grid polygons and years to process
polID.l <- c("km1", "km5", "km10", "km20")
coord.l <- c("xy1", "xy5", "xy10", "xy20")
sumby.l <- c("2007", "2008")

# Apply summaryBy function with each grid sample size for each sample year,
collapsing the lists for each year
# so that final result is a list of 3 dataframes--one for each of the five
grid sample sizes.
averageVIs <- function(setname, nres, pvec) {
# Processes VI values by year and grid cell
# Args:
#   setname: String indicating portion of name to sample
#   nres: Number of resolution to test
#   pvec: vector of column names in setname to average
# Returns:
#   Outlists of spatial points with values averaged first by year within a
cell, then across years

  # Function for counting how many years of observations are represented in
each grid cell
  ctUnique <- function(x) length(unique(x))
  formmain <- paste(paste(pvec, sep = "", collapse = " + "), "~")

  # First average within years
  out.l1 <- lapply(1:nres, function(x) {
    # summaryBy formulas for EVI and yield values by grid ID
    form <- paste(formmain, polID.l[x])  # Main averaging formula
    form2 <- paste("year ~", polID.l[x])  # Year formula

    # Process variables for each year in list
    out.sum.minor <- lapply(1:length(sumby.l), function(y) {
      print(paste("Resolution:", polID.l[x], "- Averaging observations
from", sumby.l[y]))
      ydat <- get(setname)
      ydat.yr <- ydat@data[ydat@data$year == sumby.l[y], ]
      # Run summaryBy function, passing formula from outer lapply function
      mean.vals <- summaryBy(formula(form), ydat.yr, FUN = mean, keep.names
= TRUE)
      mean.vals$year <- sumby.l[y]  # Add new column telling which year
each observation was collected
      return(mean.vals)
    })

    # Collapse list into dataframe and round to 0
    out.frame <- do.call("rbind", out.sum.minor)  # Collapse
    cind <- grep("depth|yield|yld|ndvi|evi", colnames(out.frame))
    out.frame[, cind] <- round(out.frame[, cind], 0)  # Round to 0

    # Apply year counting function in summaryBy
    yr.count <- summaryBy(formula(form2), out.frame, FUN = c(ctUnique),
keep.names = TRUE)
    out.frame$yrs4 <- ifelse(out.frame[, polID.l[x]] %in%
yr.count[yr.count$year > 3, 1], 1, 0)
    out.frame$yrs3 <- ifelse(out.frame[, polID.l[x]] %in%
yr.count[yr.count$year > 2, 1], 1, 0)
    out.frame$yrs2 <- ifelse(out.frame[, polID.l[x]] %in%
yr.count[yr.count$year > 1, 1], 1, 0)
    out.frame$yrs1 <- ifelse(out.frame[, polID.l[x]] %in%
yr.count[yr.count$year >= 1, 1], 1, 0)
    return(out.frame)
  })

  # Now average over years within polygons, according to minimum number of
years sampled in each cell
  yrs.ct.l <- list("yrs2", "yrs1")  # Year list

  # Run summaryBy to find multi-year means for NDVI metrics in each
dataframe in list over.out.sum, calculating
  # means for each possible minimum number of years of observation criterion
  out.l2 <- lapply(1:nres, function(q) {
    print(paste("Processing for resolution", polID.l[q]))
    indat <- out.l1[[q]]  # Pull dataframe corresponding to particular grid
size
    form <- paste(formmain, polID.l[q])  # Main averaging formula

    # Run function for each minimum year criterion for the dataframe
    means.by.yrs <- lapply(yrs.ct.l, function(z) {
      print(paste("Averaging for", z))
      indat.2 <- indat[indat[, as.character(z)] == 1, ]  # Subset data by n
yrs
      mean.vals2 <- summaryBy(formula(form), indat.2, FUN = mean,
keep.names = TRUE)  # Average
      cind <- grep("depth|yield|yld|ndvi|evi", colnames(mean.vals2))  # Col
indices
      mean.vals2[, cind] <- round(mean.vals2[, cind], 0)  # Round DSSAT and
VI values to 0
      mean.vals2.xy <- merge(get(coord.l[q]), mean.vals2, by.x =
polID.l[q], by.y = polID.l[q])  # Get coords
      coordinates(mean.vals2.xy) = ~x + y
      mean.vals2.xy@proj4string <- p@proj4string
      return(mean.vals2.xy)
    })
    names(means.by.yrs) <- unlist(yrs.ct.l)
    return(means.by.yrs)
  })
  names(out.l2) <- polID.l  # Name "outer" list items
  return(out.l2)
}

# String of values to average
preds <- qw(ndvi, wh.yld, wh.cv, pre, tmx, tmin, dul)
ndvi.l <- averageVIs("p", 4, preds)



On Mon, May 28, 2012 at 5:00 PM, Gavin Simpson <gavin.simpson@ucl.ac.uk>wrote:

> Dear List,
>
> I have a SpatialGrid object - a regular 5km by 5km over the entire land
> mass of Great Britain & NI - and a SpatialPointsDataFrame containing a
> @data slot with several attributes.
>
> I want to identify which points are in each grid cell and average the
> values of the individual attributes from the points within each cell.
> Not all cells will contain a point (for those an `NA` returned value is
> fine) though some cells will contain many points. Later I may wish to
> calculate some other summary (min or max) rather than mean.
>
> If I understand correctly, this gives the desired result:
>
> aggregate(fabResults.sp, by = frameGrid, FUN = mean)
>
> where `fabResults.sp` is my SpatialPointsDataFrame and `frameGrid` is my
> Spatial Grid. This takes a while to compute however.
>
> I was unable to come up with an `over()` call that worked.
>
> Have I got the `aggregate()` step right? Is there a way to use `over()`
> to do the same operation? Am I missing a better alternative?
>
> Thanks in advance,
>
> Gavin
> --
> %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
>  Dr. Gavin Simpson             [t] +44 (0)20 7679 0522
>  ECRC, UCL Geography,          [f] +44 (0)20 7679 0565
>  Pearson Building,             [e] gavin.simpsonATNOSPAMucl.ac.uk
>  Gower Street, London          [w] http://www.ucl.ac.uk/~ucfagls/
>  UK. WC1E 6BT.                 [w] http://www.freshwaters.org.uk
> %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

	[[alternative HTML version deleted]]

