检查点是否在由多个多边形/孔组成的空间对象中 [英] Check if point is in spatial object which consists of multiple polygons/holes

问题描述

我有一个带有11589个多边形"类对象的SpatialPolygonsDataFrame.这些对象中的10699个正好由1个多边形组成，但是这些对象的其余部分则由多个多边形(2至22个)组成.

I have a SpatialPolygonsDataFrame with 11589 objects of class "polygons". 10699 of those objects consists of exactly 1 polygon, however the rest of those objects consists of multiple polygons (2 to 22).

如果一个对象包含多个多边形，则可能出现三种情况:

If an object of consists of multiple polygons, three scenarios are possible:

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

1. Sometimes, those additional polygons describe a "hole" in the geographic ara describe by the first polygon in the object of class "polygons".
2. Sometimes, those additional polygons describe additional geographic areas, i.e. the shape of the region is quite complex and described by putting together multiple parts.
3. Sometimes, it might be a mix of both, 1) and 2).

Stackoverflow帮助我正确地绘制了这样的空间对象(绘制已定义的空间区域多个多边形).

Stackoverflow helped me to plot such an spatial object properly (Plot spatial area defined by multiple polygons).

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

However, I am still not able to answer how to determine whether a point (defined by longitude/latitude) is in a polygon.

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

Below is my code. I tried to apply the function point.in.polygon in the sp package, but found no way how it could handle such an object which consists of multiple polygons/holes.

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

推荐答案

您可以使用rgeos软件包中的gContains(...)轻松完成此操作.

You can do this simply with gContains(...) in the rgeos package.

gContains(sp1,sp2)

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

returns a logical depending on whether sp2 is contained within sp1. The only nuance is that sp2 has to be a SpatialPoints object, and it has to have the same projection as sp1. To do that, you would do something like this:

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

这是一个基于您上一个问题的答案的有效示例.

Here is a working example based on the answer to your previous question.

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(...)

Here, point A is in the main polygon, point B is in a lake (hole), point C is on an island, and point D is completely outside the region. So this code checks all of the points using 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

这篇关于检查点是否在由多个多边形/孔组成的空间对象中的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！