1

I am new to R, and trying to build a simple function that calculates an average (mean) of Sum_MVKT extracted from a set of predefined point locations.

For this, I have two shapefiles named Traffic_Val (polygon), and StnLocations (points).

"Traffic_Val" contains a large number of polygons and each polygon has an unique value (Sum_MVKT).

"StnLocations" consists of multiple points that lie on top of polygons in "Traffic_Val".

Provided the predefined StnLocations, I'd like to extract the value from Traffic_Val at each station location and calculate the average.

Here is the R script that I have written so far..

library(sp)
library(maptools)

setwd("C:/Users/calculateTraffic") 

# load data:
SA <- readShapePoly('Traffic_Val.shp') # This file contains a list of traffic values named as "Sum_MVKT"
proj4string(SA) <- CRS("+init=epsg:26915") 


### Existing Station Locations ###
PtsLocation <- readShapePoints('stnLocations.shp')
names(PtsLocation@data)
names(PtsLocation@data)<-c("ID", "x", "y", "Site_desc")

plot(SA)
points(PtsLocation, pch = 1, col = 3)

###########################
#### calculateTraffic ##### 
###########################

calculateTraffic = function (PtsLocation)
{
  obs = PtsLocation
  if (missing(formulaString) || is.null(formulaString)) {
    eq = dum ~ 1
  } else eq = formulaString
  ############ changed following 'if statement' 
  if (!"data.frame" %in% getSlots(class(obs)) & (terms(eq)[[3]] ==  
                                               1 || all(all.vars(eq)[-1] %in% dimnames(coordinates(obs))[[2]]))) {
    obs = SpatialPointsDataFrame(obs, data = data.frame(dum = rep(1,
                                                              dim(coordinates(obs))[1])))
    names(obs) = as.character(eq[[2]])
  }

  meanTraffic <<- mean( )

  return(meanTraffic)

}
1

Your description is a bit confusing. It sounds like you would like to aggregate the mean of values, located at points, based on the intersection with a polygon object.

Even if this is not exactly the case, this code should at least point you in the right direction.

require(sp)

# Create example data
data(meuse)
coordinates(meuse) = ~x+y
sr <- list(sr1=Polygons(list(Polygon(cbind(c(180114, 180553, 181127, 181477, 181294, 181007, 180409, 
  180162, 180114), c(332349, 332057, 332342, 333250, 333558, 333676, 
  332618, 332413, 332349)))),"1"),
sr2=Polygons(list(Polygon(cbind(c(180042, 180545, 180553, 180314, 179955, 179142, 179437, 
  179524, 179979, 180042), c(332373, 332026, 331426, 330889, 330683, 
  331133, 331623, 332152, 332357, 332373)))),"2"),
sr3=Polygons(list(Polygon(cbind(c(179110, 179907, 180433, 180712, 180752, 180329, 179875, 
  179668, 179572, 179269, 178879, 178600, 178544, 179046, 179110),
  c(331086, 330620, 330494, 330265, 330075, 330233, 330336, 330004, 
  329783, 329665, 329720, 329933, 330478, 331062, 331086)))),"3"),
sr4=Polygons(list(Polygon(cbind(c(180304, 180403,179632,179420,180304),
  c(332791, 333204, 333635, 333058, 332791)))),"4"))
srdf=SpatialPolygonsDataFrame(SpatialPolygons(sr), data.frame(row.names=c("1","2","3","4"), PIDS=1:4))

plot(srdf)
  plot(meuse, pch=20, add=TRUE)

# Use over to add polygon ID's to meuse points
meuse@data <- data.frame(meuse@data, PID=(meuse %over% srdf))

# Use tapply to calculate mean for each polygon
tapply(meuse@data$copper, meuse@data$PID, FUN=mean)
  • Thanks!!! I was able to figure it out from your example! – Tae J. Kwon Jan 23 '15 at 6:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.