28

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    <NA>
1 E05000448                   Lewisham Central    <NA>
2 E05003149                            Hawcoat    <NA>

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?

loki
  • 514
  • 4
  • 18
ben_aaron
  • 383
  • 1
  • 3
  • 6
  • Have you looked at point.in.polygon() in package sp? –  Mar 03 '15 at 22:53
  • @arvi1000 I have and will try this again. My thought about point.in.polygon was whether this would preserve the variables month and crime_type . Do you know about that? – ben_aaron Mar 03 '15 at 23:23
  • I've tried a bit more with point.in.poly and have finally selected those points that fall into the relevant polygons. Thanks. – ben_aaron Mar 04 '15 at 14:18
  • Then perhaps you should answer your own question with your solution. Remember, good answers are what this site is all about. – SlowLearner Mar 04 '15 at 14:30

3 Answers3

32

over() from package sp can be a little confusing but works well. I'm assuming you've already made "A" spatial with coordinates(A) <- ~longitude+latitude:

# Overlay points and extract just the code column: 
a.data <- over(A, B[,"code"])

Instead of a point spatial object, this simply gives you a data frame, with the same no. rows as A, and a single variable "code" from each intersecting polygon from B.

# Add that data back to A:
A$bcode <- a.data$code
Simbamangu
  • 14,773
  • 6
  • 59
  • 93
  • I have found over() to have issues with points at the vertices of the polygons, although I think this is the easiest solution I have found so far. – JMT2080AD Jul 20 '16 at 19:51
  • What issues have you had? – Simbamangu Jul 21 '16 at 10:48
  • Exclusion. I need to explore it further. I'll pm you some data later today and we can look at it together if your interested. I might be wrong, but I'm pretty sure there are some degeneracies in the algorithm that need to be taken care of, at least for my data. – JMT2080AD Jul 21 '16 at 17:33
  • Nevermind. It must be something with my data. This experimental set works fine. http://www.r-fiddle.org/#/fiddle?id=m5sTjE4N&version=1 – JMT2080AD Jul 21 '16 at 22:59
  • 2
    This is a much more straightforward approach than the accepted answer, and does not require installing additional packages other than rgdal. – Bryce Frank Feb 09 '18 at 20:35
  • @BryceFrank, actuall this is no longer the case. The spatialEco::point.in.polygon function now deals with some multipart geometries that sp::over does not handle correctly. There are argument options on exactly how results are returned from multipart overlays. – Jeffrey Evans Jun 28 '18 at 15:45
9

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)
Jeffrey Evans
  • 31,750
  • 2
  • 47
  • 94
  • @ arvi1000, yes but sp::point.in.polygon produces a logical. The spatialEco:point.in.poly is a wrapper for over but returns an sp SpatialPointsDataFrame and shortcuts some steps in relating the polygon attributes, much like raster:intersect does for rgeos::gIntersect. – Jeffrey Evans Apr 14 '16 at 16:35
  • sp::point.in.polygon actually returns a numeric value (0=point is outside, 1=inside, 2=on edge, 3=on vertex). Could be the right thing for some circumstances. Thought it was helpful to note here, since this is a top google result for related terms – arvi1000 Apr 14 '16 at 18:02
1

Here is a dplyr like solution:

library(spdplyr)

ukcounties <- geojsonio::geojson_read("data/Westminster_Parliamentary_Constituencies_December_2018_UK_BGC/uk_country.geojson",
                                      what = "sp")
pop <- read_excel("data/SAPE20DT7-mid-2017-parlicon-syoa-estimates-unformatted.xls",sheet = "data")
pop <- janitor::clean_names(pop)

ukcounties_pop <- ukcounties %>% inner_join(pop, by = c("pcon18nm" = "pcon11nm"))

The population data comes from: https://www.ons.gov.uk/peoplepopulationandcommunity/populationandmigration/populationestimates/datasets/parliamentaryconstituencymidyearpopulationestimates

I had to convert the shape files downloaded from to geoJson: https://geoportal.statistics.gov.uk/datasets/westminster-parliamentary-constituencies-december-2018-uk-bgc/data?page=1

You can do so by:

uk_constituencies <- readOGR("data/Westminster_Parliamentary_Constituencies_December_2018_UK_BGC/Westminster_Parliamentary_Constituencies_December_2018_UK_BGC.shp")
uk_constituencies # this is in tmerc format. we need to convert it to WGS84 required by geoJson format.

# First Convert to Longitude / Latitude with WGS84 Coordinate System
wgs84 = '+proj=longlat +datum=WGS84'
uk_constituencies_trans <- spTransform(uk_constituencies, CRS(wgs84))

# Convert from Spatial Dataframe to GeoJSON
uk_constituencies_json <- geojson_json(uk_constituencies_trans)

# Save as GeoJSON file on the file system.
geojson_write(uk_constituencies_json, file = "data/Westminster_Parliamentary_Constituencies_December_2018_UK_BGC/uk_country.geojson")

#read back in:
ukcounties <- geojsonio::geojson_read("data/Westminster_Parliamentary_Constituencies_December_2018_UK_BGC/uk_country.geojson",
                                      what = "sp")
Shery
  • 111
  • 2