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

Question

我有一个包含11589个类为"polygons"的对象的SpatialPolygonsDataFrame。其中10699个对象仅由一个多边形组成，但其余的对象由多个多边形（2到22个）组成。

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

1. 有时，这些附加多边形描述了类别“polygons”中第一个多边形所描述地理区域中的一个“洞”。
2. 有时，这些附加多边形描述了其他地理区域，即该地区的形状相当复杂，并通过结合多个部分来描述。
3. 有时，可能是1和2的混合。

Stackoverflow帮助我正确地绘制此类空间对象（Plot spatial area defined by multiple polygons）。

然而，我仍然不知道如何确定一个点（由经度/纬度定义）是否在一个多边形内。

以下是我的代码。 我尝试应用sp包中的point.in.polygon函数，但找不到如何处理由多个多边形/洞组成的对象的方法。

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

- majom

Answer 1

您可以在rgeos软件包中使用gContains(...)来轻松实现此操作。

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

- jlhoward

Answer 2

由于您可以使用“点在多边形内”算法，而这显然还没有被适当地设计用于处理R中的多多边形情况（我实际上觉得有点奇怪），因此您必须循环遍历每个多边形。现在的诀窍是，如果您在奇数个多边形内部，则在多重多边形内部。如果您在偶数个多边形内部，则实际上在形状外面。

使用射线交叉的点与多边形测试应该已经能够处理这个问题，只需确保将所有顶点传递给原始的点.in.polygon测试即可，但我不确定R正在使用哪种机制，因此我只能向您提供上述偶数/奇数建议。

我还找到了这段代码，不确定它是否有帮助：

require(sp)
require(rgdal)
require(maps)

# read in bear data, and turn it into a SpatialPointsDataFrame
bears <- read.csv("bear-sightings.csv")
coordinates(bears) <- c("longitude", "latitude")

# read in National Parks polygons
parks <- readOGR(".", "10m_us_parks_area")

# tell R that bear coordinates are in the same lat/lon reference system
# as the parks data -- BUT ONLY BECAUSE WE KNOW THIS IS THE CASE!
proj4string(bears) <- proj4string(parks)

# combine is.na() with over() to do the containment test; note that we
# need to "demote" parks to a SpatialPolygons object first
inside.park <- !is.na(over(bears, as(parks, "SpatialPolygons")))

# what fraction of sightings were inside a park?
mean(inside.park)
## [1] 0.1720648

# use 'over' again, this time with parks as a SpatialPolygonsDataFrame
# object, to determine which park (if any) contains each sighting, and
# store the park name as an attribute of the bears data
bears$park <- over(bears, parks)$Unit_Name

# draw a map big enough to encompass all points (but don't actually plot
# the points yet), then add in park boundaries superimposed upon a map
# of the United States
plot(coordinates(bears), type="n")
map("world", region="usa", add=TRUE)
plot(parks, border="green", add=TRUE)
legend("topright", cex=0.85,
    c("Bear in park", "Bear not in park", "Park boundary"),
    pch=c(16, 1, NA), lty=c(NA, NA, 1),
    col=c("red", "grey", "green"), bty="n")
title(expression(paste(italic("Ursus arctos"),
    " sightings with respect to national parks")))

# now plot bear points with separate colors inside and outside of parks
points(bears[!inside.park, ], pch=1, col="gray")
points(bears[inside.park, ], pch=16, col="red")

# write the augmented bears dataset to CSV
write.csv(bears, "bears-by-park.csv", row.names=FALSE)

# ...or create a shapefile from the points
writeOGR(bears, ".", "bears-by-park", driver="ESRI Shapefile")

- Ted

---

感谢您的建议。如果 point.in.polygon 能够处理多个多边形，您知道在这种情况下如何为 point.in.polygon(point.x, point.y, pol.x, pol.y, mode.checked=FALSE) 提供必要的参数吗？ - majom

我添加了一些示例代码，看起来它们使用over()函数。 - Ted

point.in.polygon 适用于简单环形多边形的底层使用。对于复杂对象的一般测试，请使用?over。 - mdsumner

谢谢！稍作修改后，我能够使用read.csv("bear-$hit-sightings.csv")和readOGR(".", "The_Woods")来回答人们多年来一直在问的另一个问题！ - geneorama