rmetrics - test if longitude and latitude coordinates are on land or sea

Question

I have a series of longitudes and latitudes for earthquakes. I would like to be able to split them into ones that are over land and ones that are over sea.

Is there an r function to do that?

Answer 1

It kinda depends on how much work you want to do.

First, grab the "oceans" shapefile from the Natural Earth site: http://www.nacis.org/naturalearth/110m/physical/ne_110m_ocean.zip

Unzip that somewhere (the example below has it on the desktop since I was being lazy), then read it in and use a function from the prevR package (which I've included below since getting prevR to install isn't worth it for just the function you need).

library(sp)
library(rgdal)
require(maptools)


# VERBATIM COPY FROM prevR package. They deserve all the credit for this function

point.in.SpatialPolygons = function(point.x, point.y, SpP){
  ###############################################################################################
  # Cette fonction renvoie pour chaque point defini par le couple (point.x, point.y) T ou F 
  #     si le point est a l'interieur ou non du spatialPolygons SpP
  # Un point est considere a l'interieur de N polygons si il est a l'interieur d'au moins
  # un polygon non Hole et a l'exterieur de tous les polygons Hole
  # Cette foncion est utilisee par toutes les fonctions de lissage (krige , kde, idw) . 
  # En effet ces fonctions travaillent sur un grid rectangulaire englobant les donnees. En presentation on ne veut que les resulats 
  #   interieurs a la frontiere qui est definie dans l'element SpP (SpP contient boundary). 
  #   Tous les elements du grid hors de la frontiere seront dans les programmes de lissage positonnes a NA
  # 
  ###############################################################################################

  X = slot(SpP,"polygons")
  is.inside = F
  for(i in 1:length(X)){
    PS   = slot(X[[i]],"Polygons")
    for(j in 1:length(PS)){
      pol.x = slot(PS[[j]],"coords")[,1]
      pol.y = slot(PS[[j]],"coords")[,2]
      pointsPosition = point.in.polygon(point.x, point.y, pol.x, pol.y)
      if(!slot(PS[[j]],"hole")) {
        is.inside = is.inside | pointsPosition != 0
      }
    }
  }
  is.outsideHole = T
  for(i in 1:length(X)){
    PS   = slot(X[[i]],"Polygons")
    for(j in 1:length(PS)){
      pol.x = slot(PS[[j]],"coords")[,1]
      pol.y = slot(PS[[j]],"coords")[,2]
      pointsPosition = point.in.polygon(point.x, point.y, pol.x, pol.y)
      if(slot(PS[[j]],"hole")) {
        is.outsideHole = is.outsideHole & (pointsPosition == 0 |  pointsPosition == 3)
      }
    }
  }
  is.inside & is.outsideHole
}


# where is the oceans' shapefile
setwd("~/Desktop/ne_110m_ocean/")

# read in the shapefile; repair=TRUE prbly isn't necessary
# but it doesn't hurt and it's a force of habit for me
oceans <- readShapePoly("ne_110m_ocean.shp", repair=TRUE)

point.in.SpatialPolygons(-105, 45, oceans)
## [1] FALSE
point.in.SpatialPolygons(-135, 30, oceans)
## [1] TRUE
point.in.SpatialPolygons(-75, -25, oceans)
## [1] TRUE
point.in.SpatialPolygons(-45, -15, oceans)
## [1] FALSE
point.in.SpatialPolygons(165, -15, oceans)
## [1] TRUE
point.in.SpatialPolygons(155, 73, oceans)
## [1] TRUE

That's not alot of testing, but you can let me know if it works for your needs.