在地图上查找多边形截距

Question

我试图在此地图上找到截取 R 中州边界的半径。

到目前为止，这是我的代码。感谢用户 Gregoire Vincke 提供大部分解决方案。

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

我想在新数据框中存储截取州边界的点，如有任何帮助，我们将不胜感激！

Answer 1

编辑：这是使用 R（sp、rgdal、rgeos）中的地理空间分析包的工作流程的广泛概述。

您不需要使用地图包和 stateMapEnv，而是需要州边界的多边形 shapefile，例如可以在此处找到的文件： https://www.census.gov/geo/maps-data/data/cbf/cbf_state.html

然后您可以使用 rgdal 包中的 readOGR 在 R 中加载该 shapefile 以获得 SpatialPolygons（我们称它为 state_poly）和一个 多边形 每个状态的对象。

根据您的 long/lat 坐标创建一个 SpatialPoints 对象：

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

此时你的pts和state_poly应该在longitude/latitude坐标，但是画圈围绕点的固定半径，您需要将它们转换为投影坐标（即以米为单位）。有关详细信息，请参阅此问题： Buffer (geo)spatial points in R with gbuffer
用每个点周围的圆圈半径创建一个矢量，并将其与 gBuffer（来自 rgeos）和点图层一起使用：

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

byid 参数意味着它对每个点单独执行，使用 radii 中的不同值，顺序与点相同。

将州多边形转换为线：state_lines <- as(state_poly, "SpatialLines")
使用gIntersects(circ, state_lines, byid = TRUE).

由于 byid = TRUE，return 值是一个矩阵，在 spgeom1 中每个圆圈一行，在 spgeom2 中每个状态边界一列。请注意，如果圆与两个状态之间的边界相交，则该行中应该有两个 "TRUE" 值（每个状态一个）。如果它与水或美国的外围相交，它可能只有一个 "TRUE" 行中的值。

Answer 2

这是最终代码！

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

maps

r

data-visualization

polygons