CHopp - 1 year ago 78

R 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!

Answer Source

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