rmetrics - 测试经度和纬度坐标是在陆地还是海上答案

【问题标题】：rmetrics - test if longitude and latitude coordinates are on land or searmetrics - 测试经度和纬度坐标是在陆地还是海上
【发布时间】：2014-06-21 00:18:18
【问题描述】：

我有一系列地震的经度和纬度。我希望能够将它们分成陆地上的和海上的。

有没有 r 函数可以做到这一点？

【问题讨论】：

标签： r gps spatial rgdal ogr

【解决方案1】：

这有点取决于你想做多少工作。

首先，从 Natural Earth 站点获取“oceans”shapefile：http://www.nacis.org/naturalearth/110m/physical/ne_110m_ocean.zip

解压到某个地方（因为我比较懒惰，下面的示例将它放在桌面上），然后将其读入并使用prevR 包中的一个函数（自从安装prevR 以来，我已将其包含在下面）仅仅为了你需要的功能是不值得的）。

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

这不是很多测试，但你可以让我知道它是否适合你的需求。

【讨论】：