CODEWITHSUNDEEP
rmapsdata-visualizationpolygons
CHopp
CHopp

Reputation: 17

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.

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!

Upvotes: 0

Views: 433

Answers (2)

CHopp
CHopp

Reputation: 17

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")

Upvotes: 0

Philippe Marchand
Philippe Marchand

Reputation: 676

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

  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.

Upvotes: 1

Related Questions