Remove points which are out of shapefile (or raster) extent

Question

My problem is quite simple. I have thousands of points (coordinates) to analyse, however some of them are so inaccurate that they "hit" the sea. I need to get rid of them. I was trying to find some R function which is capable to filter (delete) such undesirable points according to borders (mask?) of wrld_simpl (shapefile) or wrld_simpl (converted to raster)? Please only R solution.

Here is a small working example:

# coordinates
library(maptools)
library(raster)
data(wrld_simpl)
test.points <- data.frame(lon = c(-5,-3,-4,19,19,19),
                          lat = c(45,45,45,50,51,52))

# plot (shape-file + points)
plot(wrld_simpl, axes = TRUE, col = "grey",
     xlim = c(-15,+45),
     ylim = c(+35,+60)); box()
points(test.points$lon, test.points$lat,
       col = "red", pch = 20, cex = 1)

# converting wrld_simpl to raster
ext <- extent(wrld_simpl)
xy <- abs(apply(as.matrix(bbox(ext)), 1, diff))
n <- 5
r <- raster(ext, ncol=xy[1]*n, nrow=xy[2]*n)
rr <-rasterize(wrld_simpl, r)

# plot (raster + points)
plot(rr, axes = TRUE, col = "grey",
     xlim = c(-15,+45),
     ylim = c(+35,+60)); box()
points(test.points$lon, test.points$lat,
       col = "red", pch = 20, cex = 1)

Answer 1

Use over and the Natural Earth oceans shapefile:

library(rgdal)
library(rgeos)

URL <- "http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/physical/ne_110m_ocean.zip"
fil <- basename(URL)
if (!file.exists(fil)) download.file(URL, fil)
fils <- unzip(fil)
oceans <- readOGR(grep("shp$", fils, value=TRUE), "ne_110m_ocean",
                  stringsAsFactors=FALSE, verbose=FALSE)

test.points <- data.frame(lon = c(-5,-3,-4,19,19,19),
                          lat = c(45,45,45,50,51,52))

coordinates(test.points) <- ~lon+lat
proj4string(test.points) <- CRS(proj4string(oceans))

over(test.points, oceans)

##   scalerank featurecla
## 1         0      Ocean
## 2         0      Ocean
## 3         0      Ocean
## 4        NA       <NA>
## 5        NA       <NA>
## 6        NA       <NA>

Wherever there are NAs those points are not in the ocean.

Answer 2

You can do that like this:

sp <- SpatialPoints(test.points, proj4string=CRS(proj4string(wrld_simpl)))
spsel <- sp[wrld_simpl]