4

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)

enter image description here

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

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.

| improve this answer | |
1

You can do that like this:

sp <- SpatialPoints(test.points, proj4string=CRS(proj4string(wrld_simpl)))
spsel <- sp[wrld_simpl]
| improve this answer | |
  • @ Robert Hijmans Can this be done with a GADM object? I keep getting an error Error in .local(obj, ...) : cannot derive coordinates from non-numeric matrix – I Del Toro Jun 22 at 14:43
  • I am using SpatialPoints and SpatialPolygonsDataFrame objects. There is no such thing as a GADM object. – Robert Hijmans Jun 22 at 16:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.