r - rmetrics - 経度と緯度の座標が陸地か海域かをテストする

Question

私は一連の地震の経度と緯度を持っています。陸と海に分けられるようにしたいです。

それを行うr関数はありますか？

Answer 1

それは、あなたがどれだけの仕事をしたいかによって異なります。

まず、Natural Earth サイトから「海」シェープファイルを取得します: 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

これは多くのテストではありませんが、必要に応じて機能するかどうかをお知らせください。