Check if point is in spatial object which consists of multiple polygons/holes

geospatial, point-in-polygon, polygon, r, spatial

Solution

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

gContains(sp1,sp2)

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

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

Problem

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: - Sometimes, those additional polygons describe a "hole" in the geographic ara describe by the first polygon in the object of class "polygons". - 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. - Sometimes, it might be a mix of both, 1) and 2). 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. 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) ```

Original source