Join spatial point data to polygons in R

I am trying to perform a spatial join between point data and polygon data.

I have data that indicate the spatial coordinates of an event in my csv file A and have another file, shapefile B, that contains the boundaries of an area as polygons.

head(A)
  month   longitude latitude lsoa_code                   crime_type
1 2014-09 -1.550626 53.59740 E01007359        Anti-social behaviour
2 2014-09 -1.550626 53.59740 E01007359                 Public order
3 2014-09 -1.865236 53.93678 E01010646        Anti-social behaviour


head(B@data)
  code      name                                  altname
0 E05004934 Longfield, New Barn and Southfleet    
1 E05000448                   Lewisham Central    
2 E05003149                            Hawcoat    

I want to join the crime data A to my shapefile B to map the crime events that happen in my area A. Unfortunately I cannot perform an attribute join based code as the code in A refers to different units than the code in B.

I've read a number of tutorials and posts but could not find an answer. I tried:

joined = over(A, B)

and overlay, but did not accomplish what I wanted.

Is there a way to do this join directly or would an intermediate transformation from A to another format be needed?

Conceptually I want to select those points of A that fall into the code areas of B (similar to "join based on spatial location in ArcGIS").

Did someone have this issue and solved it?

Answer

The point.in.poly function in the spatialEco package returns a SpatialPointsDataFrame object of the points that intersect an sp polygon object and optionally adds the polygon attributes.

First lets add the require packages and create some example data.

require(spatialEco)
require(sp)

data(meuse)
coordinates(meuse) = ~x+y
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')
sr=SpatialPolygons(list(sr1,sr2,sr3,sr4))
srdf=SpatialPolygonsDataFrame(sr, data.frame(row.names=c('1','2','3','4'), PIDS=1:4, y=runif(4)))

Now, lets take a quick look at the data and plot it.

head(srdf@data)  # polygons
head(meuse@data) # points

plot(srdf)
points(meuse, pch=20)

Finally, we can intersect the points with the polygons. The results will be a SpatialPointsDataFrame object with, in this case, two extra attributes (PIDS, y) that were contained in the srdf polygon data.

  pts.poly <- point.in.poly(meuse, srdf)
    head(pts.poly@data)

If there is not a unique identification column in the polygon data you could easily add one.

srdf@data$poly.ids <- 1:nrow(srdf) 

Once we have the points and polygons intersected, we can aggregate the points using the unique polygon ID's that were an attribute in the polygon data.

# Number of points in each polygon
tapply(pts.poly@data$lead, pts.poly@data$PIDS, FUN=length)

# Mean lead in each polygon
tapply(pts.poly@data$lead, pts.poly@data$PIDS, FUN=mean)