Така че тук се случват няколко неща. Първо, вашият набор от данни изглежда има рН спрямо дълбочина. Така че докато има ~ 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. След това трансформираме и двете в проекцията на 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))
РЕДАКТИРАНЕ: Коментарът на 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
