Find Polygon Intercepts on a Map

Question

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

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

library("maps")
library("mapproj")
library("RColorBrewer")
library("mapdata")
library("ggplot2")
library("rgeos")
library("dismo")
library("ggmap")
library("rgdal")

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

plotCircle(-71.4868,42.990684,20)
plotCircle(-72.57085,41.707932,12)
...
#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!

score 1 · Answer 1 · edited May 23 '17 at 12:01

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

Instead of using the maps package and stateMapEnv, you want a polygon shapefile of state boundaries, like one that can be found here: https://www.census.gov/geo/maps-data/data/cbf/cbf_state.html

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.

Create a SpatialPoints object from your long/lat coordinates:

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

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
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.

Convert the state polygons to lines: state_lines <- as(state_poly, "SpatialLines")
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.

Although intersection on a water or cnadian boundary would not meet your test. I was thinking you might be able to compare the area of the circle with the area of the intersection and if different that implies a boundary crossing. — IRTFM, Aug 03 '16 at 21:30
That's true. You could also convert the states to SpatialLines in which case the circles would not intersect if fully inside boundaries. — Philippe Marchand, Aug 04 '16 at 01:28
Hi Phillippe, thanks for the advice, just not sure how to use the arguments to check all polygons? ie spgeom1 and 2 — CHopp, Aug 04 '16 at 13:14
@CHopp All the polygons? You didn't offer a method to make _any_ polygons. — IRTFM, Aug 04 '16 at 15:19
I modified my answer above to list all the steps, which includes creating the spatial objects and projecting them (because you can't easily do fixed-radius circles in lat/long coordinates). — Philippe Marchand, Aug 04 '16 at 17:12

CHopp · Accepted Answer · 2016-08-15T16:06:31.397

Here is the Final Code!

library("maps")
library("mapproj")
library("RColorBrewer")
library("mapdata")
library("ggplot2")
library("rgeos")
library("dismo")
library("ggmap")
library("rgdal")

#import shape file (.shp), make sure all the other files in the zip are included in
#your file location!
state_poly <- readOGR(dsn = 'C:/Users/chopp/Documents/R', layer='cb_2015_us_state_500k')

#data containing lng and lat coordinates with radii
data <- read.csv("R/longlat.csv", header = T)

#create spatial point objects out of your lng and lat data
pts <- SpatialPoints(data[,c("lng","lat")], proj4string = CRS("+proj=longlat"))

#convert spatial points to projected coordinates (points and map lines)
ptsproj <- spTransform(pts, CRS("+init=epsg:3347"))
state_poly_proj<- spTransform(state_poly, CRS("+init=epsg:3347"))

#convert radii units to meters, used in our gBuffer argument later on 
radii <- data$rad*1609.344

#create circular polygons with. byid = TRUE will create a circle for each point
circ <- gBuffer(ptsproj, width = radii, byid = TRUE)



#convert state polygons to state lines    
state_lines<- as(state_poly_proj, "SpatialLines")

#use gIntersects with byid = TRUE to return a matrix where "TRUE" represents 
#crossing state boundaries or water
intdata <- gIntersects(circ, state_lines, byid = TRUE)

#write the matrix out into a csv file
write.csv(intdata,"R/Agents Intercepts 2.csv")

Find Polygon Intercepts on a Map

2 Answers2