океан широта долгота точка расстояние от берега

Question

Я начал «свободный» проект с открытым исходным кодом, чтобы создать новый набор данных для рН земных океанов.

Я начал с открытой набор данных из NOAA и создал 2,45 миллионов строк набора данных с теми столбцами:

colnames(NOAA_NODC_OSD_SUR_pH_7to9) 
[1] "Year" "Month" "Day" "Hour" "Lat" "Long" "Depth" "pH" 

Метод документа HERE.

Data-set HERE.

Моя цель теперь состоит в том, чтобы «квалифицировать» каждую строку (2.45 м) ... для этого мне нужно рассчитать расстояние от каждой точки Lat/Long до ближайшего берега.

Так что я ищу для метода, который будет принимать В: Lat/Long Out: Расстояние (км от берега)

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

У меня есть поиск метода для этого, но все, кажется, нужны пакеты/программное обеспечение, которых у меня нет.

Если кто-то захочет помочь, я был бы признателен. Или, если вы знаете простой (бесплатно) метод для достижения этой цели, пожалуйста, дайте мне знать ...

Я могу работать в R программирование, Shell скрипты вещи, но не эксперт из тех ....

Answer 1

Итак, здесь есть несколько вещей. Во-первых, ваш набор данных, по-видимому, имеет значение pH и глубину. Так что, когда есть строки ~ 2.5MM, есть только ~ 200 000 строк с глубиной = 0 - все еще много.

Во-вторых, чтобы добраться до ближайшего побережья, вам нужен шейп-файл береговых линий. К счастью, это доступно here, на отличном Natural Earth website.

В-третьих, ваши данные длинные/лат (так, единицы = градусы), но вы хотите расстояние в км, поэтому вам нужно преобразовать свои данные (данные береговой линии выше также в long/lat, а также быть трансформированным). Одна из проблем с преобразованиями состоит в том, что ваши данные, очевидно, глобальны, и любое глобальное преобразование обязательно будет непланарным. Таким образом, точность будет зависеть от фактического местоположения. Правильный способ сделать это состоит в том, чтобы скомпоновать ваши данные, а затем использовать набор планарных преобразований, соответствующих любой сетке, в которой находятся ваши точки. Однако это выходит за рамки этого вопроса, поэтому мы будем использовать глобальное преобразование (mollweide) просто чтобы дать вам представление о том, как это делается в R.

library(rgdal) # for readOGR(...); loads package sp as well 
library(rgeos) # for gDistance(...) 

setwd(" < directory with all your files > ") 
# WGS84 long/lat 
wgs.84 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" 
# ESRI:54009 world mollweide projection, units = meters 
# see http://www.spatialreference.org/ref/esri/54009/ 
mollweide <- "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" 
df  <- read.csv("OSD_All.csv") 
sp.points <- SpatialPoints(df[df$Depth==0,c("Long","Lat")], proj4string=CRS(wgs.84)) 

coast <- readOGR(dsn=".",layer="ne_10m_coastline",p4s=wgs.84) 
coast.moll <- spTransform(coast,CRS(mollweide)) 
point.moll <- spTransform(sp.points,CRS(mollweide)) 

set.seed(1) # for reproducible example 
test <- sample(1:length(sp.points),10) # random sample of ten points 
result <- sapply(test,function(i)gDistance(point.moll[i],coast.moll)) 
result/1000 # distance in km 
# [1] 0.2185196 5.7132447 0.5302977 28.3381043 243.5410571 169.8712255 0.4182755 57.1516195 266.0498881 360.6789699 

plot(coast) 
points(sp.points[test],pch=20,col="red") 

Так что читает набор данных, извлекает строки, где Depth==0, и преобразует его в объект через SpatialPoints. Затем мы читаем базу данных береговых линий, загруженную со ссылкой выше, в объект SpatialLines. Затем мы преобразуем как проекцию Mollweide с использованием spTransform(...), затем используем gDistance(...) в пакете rgeos для расчета минимального расстояния между каждой точкой и ближайшим побережьем.

Опять же, важно помнить, что, несмотря на все десятичные знаки, эти расстояния приблизительные.

Одна очень большая проблема - это скорость: этот процесс занимает ~ 2 минуты на 1000 расстояний (в моей системе), поэтому для запуска всего 200 000 расстояний потребуется около 6,7 часов. Одним из вариантов, теоретически, было бы найти базу данных береговой линии с более низким разрешением.

Код ниже рассчитает все 201 000 дистанций.

## not run 
## estimated run time ~ 7 hours 
result <- sapply(1:length(sp.points), function(i)gDistance(sp.points[i],coast)) 

EDIT: OP Комментарий о сердечников заставил меня думать, что это может быть случай, когда улучшение от распараллеливания может быть стоит усилий. Итак, вот как вы могли бы запустить это (в Windows) с помощью параллельной обработки.

library(foreach) # for foreach(...) 
library(snow)  # for makeCluster(...) 
library(doSNOW) # for resisterDoSNOW(...) 

cl <- makeCluster(4,type="SOCK") # create a 4-processor cluster 
registerDoSNOW(cl)    # register the cluster 

get.dist.parallel <- function(n) { 
    foreach(i=1:n, .combine=c, .packages="rgeos", .inorder=TRUE, 
      .export=c("point.moll","coast.moll")) %dopar% gDistance(point.moll[i],coast.moll) 
} 
get.dist.seq <- function(n) sapply(1:n,function(i)gDistance(point.moll[i],coast.moll)) 

identical(get.dist.seq(10),get.dist.parallel(10)) # same result? 
# [1] TRUE 
library(microbenchmark) # run "benchmark" 
microbenchmark(get.dist.seq(1000),get.dist.parallel(1000),times=1) 
# Unit: seconds 
#      expr  min  lq  mean median  uq  max neval 
#  get.dist.seq(1000) 140.19895 140.19895 140.19895 140.19895 140.19895 140.19895  1 
# get.dist.parallel(1000) 50.71218 50.71218 50.71218 50.71218 50.71218 50.71218  1 

Использование 4 ядер повышает скорость обработки примерно в 3 раза Так, с 1000 расстояний занимает около минуты, 100000 должны занять немного меньше, чем за 2 часа.

Обратите внимание, что использование times=1 является злоупотреблением microbenchmark(...) действительно, так как все дело в том, чтобы запустить процесс несколько раз и усреднить результаты, но у меня просто не было терпения.