r - 检查点是否在由多个多边形/孔组成的空间对象中

Question

我有一个带有11589类“多边形”的对象的SpatialPolygonsDataFrame。这些对象中的10699个正好由1个多边形组成，但是其余这些对象则由多个多边形（2至22个）组成。

如果一个对象由多个多边形组成，则可能出现以下三种情况：


有时，这些其他多边形在地理区域中描述了一个“孔”，该地理区域由“多边形”类的对象中的第一个多边形描述。
有时，这些其他多边形描述了其他地理区域，即该区域的形状非常复杂，并且通过将多个部分放在一起来描述。
有时，它可能是1）和2）的混合。


Stackoverflow帮助我正确地绘制了这样的空间对象（Plot spatial area defined by multiple polygons）。

但是，我仍然无法回答如何确定点（由经度/纬度定义）是否在多边形中。

下面是我的代码。我尝试在point.in.polygon包中应用函数sp，但没有找到如何处理包含多个多边形/孔的对象的方法。

# Load packages
# ---------------------------------------------------------------------------
library(maptools)
library(rgdal)
library(rgeos)
library(ggplot2)
library(sp) 


# Get data
# ---------------------------------------------------------------------------
# Download shape information from the internet
URL <- "http://www.geodatenzentrum.de/auftrag1/archiv/vektor/vg250_ebenen/2012/vg250_2012-01-01.utm32s.shape.ebenen.zip"
td <- tempdir()
setwd(td)
temp <- tempfile(fileext = ".zip")
download.file(URL, temp)
unzip(temp)

# Get shape file
shp <- file.path(tempdir(),"vg250_0101.utm32s.shape.ebenen/vg250_ebenen/vg250_gem.shp")

# Read in shape file
map <- readShapeSpatial(shp, proj4string = CRS("+init=epsg:25832"))

# Transform the geocoding from UTM to Longitude/Latitude
map <- spTransform(map, CRS("+proj=longlat +datum=WGS84"))


# Pick an geographic area which consists of multiple polygons
# ---------------------------------------------------------------------------
# Output a frequency table of areas with N polygons 
nPolys <- sapply(map@polygons, function(x)length(x@Polygons))

# Get geographic area with the most polygons
polygon.with.max.polygons <- which(nPolys==max(nPolys))

# Get shape for the geographic area with the most polygons
Poly.coords <- map[which(nPolys==max(nPolys)),]


# Plot
# ---------------------------------------------------------------------------
# Plot region without Google maps (ggplot2) 
plot(Poly.coords,  col="lightgreen")


# Find if a point is in a polygon 
# ---------------------------------------------------------------------------
# Define points 
points_of_interest <- data.frame(long=c(10.5,10.51,10.15,10.4), 
                     lat =c(51.85,51.72,51.81,51.7),
                     id  =c("A","B","C","D"), stringsAsFactors=F)

# Plot points
points(points_of_interest$long, points_of_interest$lat, pch=19)

Answer 1

您可以使用gContains(...)包中的rgeos轻松执行此操作。

gContains(sp1,sp2)

根据sp2是否包含在sp1中，返回一个逻辑。唯一的区别是sp2必须是SpatialPoints对象，并且必须具有与sp1相同的投影。为此，您需要执行以下操作：

point <- data.frame(lon=10.2, lat=51.7)
sp2   <- SpatialPoints(point,proj4string=CRS(proj4string(sp1)))
gContains(sp1,sp2)

这是一个基于上一个问题的答案的工作示例。

library(rgdal)   # for readOGR(...)
library(rgeos)   # for gContains(...)
library(ggplot2)

setwd("< directory with all your files >")
map <- readOGR(dsn=".", layer="vg250_gem", p4s="+init=epsg:25832")
map <- spTransform(map, CRS("+proj=longlat +datum=WGS84"))

nPolys <- sapply(map@polygons, function(x)length(x@Polygons))
region <- map[which(nPolys==max(nPolys)),]

region.df <- fortify(region)
points <- data.frame(long=c(10.5,10.51,10.15,10.4), 
                     lat =c(51.85,51.72,51.81,51.7),
                     id  =c("A","B","C","D"), stringsAsFactors=F)

ggplot(region.df, aes(x=long,y=lat,group=group))+
  geom_polygon(fill="lightgreen")+
  geom_path(colour="grey50")+
  geom_point(data=points,aes(x=long,y=lat,group=NULL, color=id), size=4)+
  coord_fixed()

在这里，点A在主多边形中，点B在湖（洞）中，点C在岛上，点D完全在区域之外。因此，此代码使用gContains(...)检查所有点

sapply(1:4,function(i)
  list(id=points[i,]$id,
       gContains(region,SpatialPoints(points[i,1:2],proj4string=CRS(proj4string(region))))))
#    [,1] [,2]  [,3] [,4] 
# id "A"  "B"   "C"  "D"  
#    TRUE FALSE TRUE FALSE