海洋纬度经度点距离岸边

Question

我开始了一个"免费"的开源项目,为地球海洋的pH值创建一个新的数据集.

我从NOAA的开放数据集开始,使用这些列创建了一个2.45百万行的数据集:

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

方法文件HERE.

数据集在这里.

我现在的目标是"限定"每一行(2.45米)......为此,我需要计算从纬度/长度的每个点到最近岸的距离.

所以我正在寻找一种方法:拉/长出:距离(距岸边公里)

有了这个,我可以确定数据点是否会受到岸上污染的影响,例如附近的城市污水.

我已经找到了一种方法来做到这一点,但似乎都需要我没有的软件包/软件.

如果有人愿意帮忙,我将不胜感激.或者,如果你知道一个简单(免费)的方法来实现这一目标,请告诉我...

我可以在R编程,Shell脚本中工作,但不是那些专家....

Answer 1

所以这里有几件事情.首先,您的数据集似乎具有pH与深度.因此,虽然有大约2.5MM的行,但只有大约200,000行,深度= 0 - 仍然很多.

其次,为了到达最近的海岸,你需要一个海岸线的形状文件.幸运的是,这里可以在优秀的自然地球网站上找到.

三,你的数据在长/ 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))

编辑: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个距离大约需要1分钟,因此100,000个应该花费不到2个小时.

请注意,使用times=1是microbenchmark(...)真的滥用,因为重点是多次运行该过程并平均结果,但我只是没有耐心.