CHopp CHopp - 2 months ago 5x
R Question

Find Polygon Intercepts on a Map

I'm trying to find the Radii on this map that intercept state borders in R.
enter image description here

Here is my code so far. Thanks to user Gregoire Vincke for providing much of the solution.


data("stateMapEnv") #US state map

dat <- read.csv("R/longlat.csv",header = T)

map('state',fill = T, col = brewer.pal(9,"Pastel2"))

#draws circles around a point, given lat, long and radius

plotCircle <- function(lonDec, latDec, mile) {
ER <- 3959
angdeg <- seq(1:360)
lat1rad <- latDec*(pi/180)
lon1rad <- lonDec*(pi/180)
angrad <- angdeg*(pi/180)
lat2rad <- asin(sin(lat1rad)*cos(mile/ER) + cos(lat1rad)*sin(mile/ER)*cos(angrad))
lon2rad <- lon1rad + atan2(sin(angrad)*sin(mile/ER)*cos(lat1rad),cos(mile/ER)-sin(lat1rad)*sin(lat2rad))
lat2deg <- lat2rad*(180/pi)
lon2deg <- lon2rad*(180/pi)
polygon(lon2deg,lat2deg,lty = 1 , col = alpha("blue",0.35))

point <- mapproject(dat$lng,dat$lat)
points(point, col = alpha("black",0.90), cex = 0.4, pch = 20) #plots points

#this goes on for every point

I want to store the points that intercept state borders in a new data frame, any help would be appreciated!


EDIT: Here's a broad overview of the workflow using the geospatial analyses packages in R (sp, rgdal, rgeos).

  1. Instead of using the maps package and stateMapEnv, you want a polygon shapefile of state boundaries, like one that can be found here:

You can then load that shapefile in R with readOGR from the rgdal package to get a SpatialPolygons (let's call it state_poly) with one Polygons object per state.

  1. Create a SpatialPoints object from your long/lat coordinates:

pts <- SpatialPoints(dat[, c("lng", "lat")], proj4string = CRS("+proj=longlat"))

  1. At this point your pts and state_poly should be in longitude/latitude coordinates, but to draw circles of a fixed radius around points, you need to convert them to projected coordinates (i.e. in meters). See this question for more details: Buffer (geo)spatial points in R with gbuffer

  2. Create a vector with the radii of your circles around each point, and use it with gBuffer (from rgeos) and your points layer:

circ <- gBuffer(pts, width = radii, byid = TRUE)

The byid argument means it does it separately for each point, using the different values in radii in the same order as the points.

  1. Convert the state polygons to lines: state_lines <- as(state_poly, "SpatialLines")

  2. Use gIntersects(circ, state_lines, byid = TRUE) .

Because of byid = TRUE, the return value is a matrix with one row per circle in your spgeom1 and one column per state boundaries in spgeom2. Note that if the circle intersect a boundary between two states, it should have two "TRUE" values in that row (one for each state). If it intersects with water or the external perimeter of the US it may have only one "TRUE" value in the row.