8

In R, how do I 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)

I would like to create a column called "Land" of type "logical".

15

Here in 2020, an answer that uses the sf package:

library(sf)
library(spData) ## For `world`, an sf MULTIPOLYGON object

## Create an sf POINTS 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 <- st_as_sf(points, coords=1:2, crs=4326)

## Find which points fall over land
ii <- !is.na(as.numeric(st_intersects(pts, world)))

## Check that it worked
plot(st_geometry(world))
plot(pts, col=1+ii, pch=16, add=TRUE)

enter image description here


Original 2012 answer, using the sp package:

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

| improve this answer | |

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.