Subset Spatial Points to Extract those inside of a Polygon (Countries Borders)

Question

I have a data.frame with lat long coordinates:

df<-data.frame(
     lat=c(40, 30, 40.864),
     lon=c(0, 20, 1.274)
    )

And the border of a country (Spain),

library(raster)
border <- getData("GADM",country="Spain",level=0)

I would like to select for df only the points inside the border.

How can I do this?

Note: in my reproducible example, df, the first entry point is inside, the second is clearly outside and the third is outside but close to the coast.

https://stackoverflow.com/questions/21971447/check-if-point-is-in-spatial-object-which-consists-of-multiple-polygons-holes — M--, Dec 14 '18 at 20:31
In fact the question @Masoud links is very similar to my question. However, it does not work. sapply(1:nrow(df), function(i){gContains(spa.mask, SpatialPoints(df[i,], proj4string = CRS(proj4string(spa.mask))))}) This returns `FALSE` for the three points. I do not understand why. Maybe because my border contains several polygons (islands and peninsula)? — Javi_VM, Dec 17 '18 at 08:39
As I say in the question, one of them is inside the border. I do not underestand what you mean by "map". I made a plot and checked that the point falls inside the polygon and the other two outside the polygon... — Javi_VM, Dec 19 '18 at 18:39
Check my answer. My guess is that you were implementing part of that solution in a wrong way. — M--, Dec 19 '18 at 19:55

score 1 · Accepted Answer · answered Dec 19 '18 at 19:54

While I believe the post that I referred to answers this question, I think some clarification is needed; hence, posting an answer.

To find the intersection of two geographic feature, we need to have the same projection. sp library helps us to do so. Make sure to put longitude and latitude in the right place:

sp::SpatialPoints(c(my_point$long,my_point$lat),proj4string=CRS(proj4string(my_raster)))

Using rgeos library we can check if there is an intersection between two spatial dataset/feature:

rgeos::gContains(my_raster,my_projected_point)

So, here is how it works for the OP's example:

library(sp)        #for projection
library(raster)    #for getting the border data
library(rgeos)     #for finding intersection
library(tidyverse) #for illustration only

#data
border <- getData("GADM",country="Spain",level=0)
df <- data.frame(
 lat=c(40, 30, 40.864),
 lon=c(0, 20, 1.274)
                )

#this is the part that actually check if a point is inside the border
#adapted from https://stackoverflow.com/questions/21971447
sapply(1:3,function(i)
  list(id=i,
       intersect= gContains(border,SpatialPoints(df[i,2:1],proj4string=CRS(proj4string(border))))))

#           [,1] [,2]  [,3] 
# id        1    2     3    
# intersect TRUE FALSE FALSE

#a map for better understanding
ggplot()+
  geom_polygon(data=border, aes(x=long, y=lat, group=group), 
                       fill=NA, color="grey50", size=0.25) +
  geom_point(data=df,aes(x=lon,y=lat), color="red", size=1)

enter image description here

Subset Spatial Points to Extract those inside of a Polygon (Countries Borders)

1 Answers1