Итак, здесь происходит несколько вещей. Во-первых, в вашем наборе данных, кажется, есть зависимость pH от глубины. Таким образом, в то время как есть ~ 2,5 мм строк, есть только ~ 200 000 строк с глубиной = 0 - все еще много.
Во-вторых, чтобы получить расстояние до ближайшего побережья, вам нужен шейп-файл береговых линий. К счастью, он доступен здесь, на отличном сайте Natural Earth.
В-третьих, ваши данные представлены в долготе/широте (таким образом, единицы измерения = градусы), но вам нужно расстояние в км, поэтому вам нужно преобразовать ваши данные (приведенные выше данные береговой линии также представлены в долготе/широте и также должны быть преобразованы). Одна из проблем с преобразованиями заключается в том, что ваши данные явно глобальны, а любое глобальное преобразование обязательно будет непланарным. Таким образом, точность будет зависеть от фактического местоположения. Правильный способ сделать это - создать сетку ваших данных, а затем использовать набор планарных преобразований, соответствующих той сетке, в которой находятся ваши точки. Однако это выходит за рамки этого вопроса, поэтому мы будем использовать глобальное преобразование (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. Затем мы преобразуем оба в проекцию Моллвейде, используя 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 расстояний занимает около минуты, 100 000 должно занимать чуть меньше 2 часов.
Обратите внимание, что использование times=1 на самом деле является злоупотреблением microbenchmark(...), поскольку весь смысл в том, чтобы запустить процесс несколько раз и усреднить результаты, но у меня просто не хватило терпения.
person
jlhoward
schedule
30.12.2014
