# Given a Lat and Lon, identify if the point is over land or ocean?

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"

## migrated from stackoverflow.comOct 21 '13 at 19:37

This question came from our site for professional and enthusiast programmers.

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

• 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