7

In R, what is the quickest way to identify if a point is over land or ocean?

e.g. for:

set.seed(0)
lat <- runif(10, -80, 80)
lon <- runif(10, -180, 180)
points <- expand.grid(lat, lon)
  • How can I create a column called "Land" of type "logical"
14

Something like this should work:

## One example of a SpatialPolygons object mapping Earth's land areas
library(maptools)
data(wrld_simpl)

## Create a SpatialPoints object
set.seed(0)
lat <- runif(10, -80, 80)
lon <- runif(10, -180, 180)
points <- expand.grid(lon, lat)  # Note that I reversed OP's ordering of lat/long
pts <- SpatialPoints(points, proj4string=CRS(proj4string(wrld_simpl)))

## Find which points fall over land
ii <- !is.na(over(pts, wrld_simpl)$FIPS)

## Check that it worked
plot(wrld_simpl)
points(pts, col=1+ii, pch=16)

enter image description here

  • Thanks! That is much easier than the wormhole I got into trying to mask netCDF files directly using nco although nco makes it look easy, it is difficult to understand. – Abe Oct 4 '13 at 5:22
  • Glad to help. Turns out colSums(rgeos::gIntersects(pts, wrld_simpl, byid=TRUE)) also works (& see too rgeos::gIntersection() for something slightly different). If you've got only a few ocean pts, you can get their distance from land with geosphere::dist2Line() which takes a SpatialPoints* and a SpatialPolygons*. It's really slow though, taking something like 20 seconds per SpatialPoint when wrld_simpl is the SpatialPolygons object. Might be faster if you reduced the # of component polygons by merging whichever you can, though prob. better to just find a different approach. – Josh O'Brien Oct 4 '13 at 16:23

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.