Проверьте, находится ли точка в пространственном объекте, который состоит из нескольких полигонов/отверстий

Ответ 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([email protected], function(x)length([email protected]))
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

Ответ 2

Поскольку вы можете использовать подпрограмму "point in polygon", и это, по-видимому, уже не подходит для обработки многополигонального футляра в R (который я нахожу немного нечетным на самом деле), вам остается выполнить цикл через каждый из многоугольников. Теперь трюк, если вы находитесь внутри нечетного числа полигонов, вы находитесь внутри многополигона. Если вы находитесь внутри четного числа полигонов, вы фактически находитесь вне формы.

Точка в полигональном тестировании, использующая перекрестку лучей, должна УЖЕ быть в состоянии справиться с этим, просто убедившись, что вы передаете все вершины в исходный тест point.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")