3

I have a short R script in which I read two sets of data:

  • a shapefile showing the zip codes of a European country

  • an Excel showing multiple points coordinates throughout the country

I used the points to define buffers around them. These buffers naturally intersect with the zip codes polygons and one polygon can intersect multiple buffers. Here is the code so far:

F_A <- read.xls("F_A.xlsx", sheet = "FA")
F_A$lon <- as.numeric(levels(F_A$lon))[F_A$lon]
F_A$lat <- as.numeric(levels(F_A$lat))[F_A$lat]
F_A <- na.omit(F_A)

coords <- cbind(F_A$lon, F_A$lat)
pts <- SpatialPointsDataFrame(coords = coords, data = F_A, proj4string = CRS("+proj=longlat +datum=WGS84"))
pts_proj <- spTransform(pts, CRS("+init=epsg:3347"))

Buffers <- gBuffer(pts_proj, width = 7000, byid = TRUE)
Buffers_xy <- spTransform(Buffers, CRS("+proj=longlat +datum=WGS84"))

zips <- readOGR(dsn = "Z", layer = "PLZ")
zips.wgs84 <- spTransform(zips, CRS("+proj=longlat +datum=WGS84"))
zips.wgs84.df <- fortify(zips.wgs84)

Polygons <- readOGR(dsn = "Poly", layer = "POLYGON")
Polygons.wgs84 <- spTransform(Polygons, CRS("+proj=longlat +datum=WGS84"))
Polygons.wgs84.df <- fortify(Polygons.wgs84)

I would like to add a new attribute to the zip code object counting the number of buffers intersecting with each zip code but fail to find the right command. How would one go about doing this?

To be complete, that is a follow up question on a somewhat similar issue I had already asked about (QGIS creating a heatmap) and managed to solve using QGIS.

  • Welcome to GIS SE. As a new user please take the Tour, which explains our policies on coding questions. GIS SE is less of a place where you tell us what you want, and more of a place where you show us what you've done, and where you're stuck. – Vince Jan 8 '17 at 15:28
  • My bad, apologies, i took the tour, thanks for pointing that out. I modified the original post. – Romain Jan 8 '17 at 16:54
2

You want gIntersects.

Here's a reproducible example:

#gIntersects is in rgeos; I'm only using maptools for the
#  shapefile for North Carolina
library(maptools)
library(rgeos)
NC <- readShapeSpatial(system.file("shapes/sids.shp", package="maptools")[1L],
 IDvar="FIPSNO", proj4string=CRS("+proj=longlat +ellps=clrk66"))
#gBuffer requires a CRS that has easy-to-understand units (this one has feet)
NC = spTransform(NC, CRS("+init=epsg:2264"))

#pick 10 points at random within the state
pts = spsample(NC, 10, type = "random")

#buffer (width chosen by trial&error)
pts_buffer = gBuffer(pts, width = 1e5, byid = TRUE)

#gIntersects will return an m x n logical matrix
#  m: length(NC); n: length(pts_buffer)
#  saying whether the polygons designated in each row/column intersect
#rowSums counts within each NC county how many such intersections occur
NC$count_buffers = rowSums(gIntersects(pts_buffer, NC, byid = TRUE))

#plot for confirmation
png("~/Desktop/NC.png")
plot(NC, main = "North Carolina with Random Points & Buffers")
plot(pts, add = TRUE, pch = "*", col = "blue")
plot(pts_buffer, add = TRUE, border = "blue")
dev.off()

enter image description here

| improve this answer | |

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.