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

У меня есть SpatialPolygonsDataFrame с 11589 объектами класса «полигоны». 10699 из этих объектов состоят ровно из 1 многоугольника, однако остальные объекты состоят из нескольких многоугольников (от 2 до 22).

Если объект состоит из нескольких полигонов, возможны три сценария:

Иногда эти дополнительные полигоны описывают «дыру» в географической области, описываемой первым полигоном в объекте класса «полигоны».Иногда эти дополнительные полигоны описывают дополнительные географические области, то есть форма региона довольно сложна и описывается путем объединения нескольких частей.Иногда это может быть сочетание 1) и 2).

Stackoverflow помог мне правильно построить такой пространственный объект (Пространственная площадь участка, определяемая несколькими полигонами).

Однако я все еще не могу ответить, как определить, находится ли точка (определенная по долготе / широте) в многоугольнике.

Ниже мой код.Я пытался применить функцию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([email protected], function(x)length([email protected]))

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

Ответы на вопрос(2)

Ваш ответ на вопрос