海洋纬度经度点距离海岸的距离

Question

海洋纬度经度点距离海岸的距离

rgoogle-mapsgeolocationlatitude-longitude

6

我开始了一个“免费”的开源项目，旨在创建地球海洋pH值的新数据集。

我从NOAA的开放数据集开始，并创建了一个包含245万行数据的数据集，其中包含以下列：

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

方法文档在这里。

数据集在这里。

我的目标是“确认”每个行（2.45m）......为此，我需要计算每个Lat / Long点到最近海岸的距离。

因此，我正在寻找一种方法，可以采取输入：Lat / Long 输出：距离（km离岸）

通过这样做，我可以确定数据点是否可能受到岸边污染的影响，例如附近的城市流出物。

我已经搜索了一种方法来完成这项工作，但似乎都需要我没有的软件包/软件。

如果有人愿意帮忙，我会感激。

或者，如果您知道一个简单（免费）的方法来完成此操作，请告诉我...

我可以使用R编程，Shell脚本等，但不是这些方面的专家....

- Simon Filiatrault

1

这个有用吗？（https://dev59.com/BYXca4cB1Zd3GeqPPOhE#27391421）还是这个？（https://dev59.com/Fnvaa4cB1Zd3GeqPFqGv#21302609） - jlhoward

好的，从这里看来，在R中有一些方法可以实现这个。我会继续阅读，但我远远没有完全理解这个。我希望有人能帮助我，但如果不可能，我可以自学！谢谢！ - Simon Filiatrault

你可以考虑在http://gis.stackexchange.com/上发布这个问题。 - jlhoward

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- jlhoward · Accepted Answer

所以这里有几件事情要处理。首先，您的数据集似乎具有pH vs.深度。因此，虽然有约2.5百万行，但只有约20万行深度为0-仍然很多。

其次，要获取到最近海岸的距离，您需要一个海岸线的shapefile。幸运的是，这可以在优秀的Natural Earth website上找到here。

第三，你的数据是以经纬度表示的（因此单位为度），但你需要以公里表示距离，所以你需要对数据进行转换（上面的海岸线数据也是以经纬度表示的，同样需要转换）。转换的一个问题在于，你的数据显然是全球性的，而任何全球性的转换都必然是非平面的。因此精度将取决于实际位置。正确的做法是对数据进行网格化处理，然后使用适用于点所在网格的一组平面变换。虽然如此，但这超出了本问题的范围，因此我们将使用一个全局变换（莫尔韦德投影），只是为了让你了解在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))

编辑: 根据楼主有关核心的评论，我开始思考这可能是一个值得并行化改进的实例。因此，以下是在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(...)的滥用，因为整个过程的重点是运行多次并平均结果，但我没有耐心等待。