所以这里发生了几件事。首先,您的数据集似乎具有 pH 与深度。因此,虽然有 ~ 2.5MM 行,但深度 = 0 的行只有 ~200,000 行 - 仍然很多。
其次,要获得到最近海岸的距离,您需要一个海岸线的 shapefile。幸运的是,这在 here 上可用,在出色的 Natural Earth website 上。
第三,你的数据是long/lat(所以,单位=度),但是你想要km的距离,所以你需要转换你的数据(上面的海岸线数据也是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 对象。然后我们使用spTransform(...)将两者都转换为Mollweide投影,然后我们使用rgeos包中的gDistance(...)来计算每个点与最近海岸之间的最小距离。
同样,重要的是要记住,尽管有所有小数位,这些距离只是近似值。
一个非常大的问题是速度:这个过程大约需要 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(...) 的滥用,因为重点是多次运行该过程并平均结果,但我只是没有耐心。