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

Это не много испытаний, но вы можете сообщить мне, если это работает для ваших нужд.