【问题标题】:ocean latitude longitude point distance from shore海洋 纬度 经度 点 离岸 距离
【发布时间】:2015-02-26 03:34:38
【问题描述】:

我启动了一个“免费”的开源项目来创建一个新的地球海洋 pH 值数据集。

我从 NOAA 的开放数据集开始,创建了一个包含这些列的 245 万行数据集:

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

方法文档HERE

数据集HERE.

我现在的目标是“限定”每一行(2.45m)...为此,我需要计算从纬度/经度的每个点到最近海岸的距离。

所以我正在寻找一种方法 在:纬度/经度 出:距离(离岸公里)

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

我正在寻找一种方法来做到这一点,但似乎都需要我没有的软件包/软件。

如果有人愿意提供帮助,我将不胜感激。 或者,如果您知道实现此目的的简单(免费)方法,请告诉我...

我可以从事 R 编程、Shell 脚本方面的工作,但不是这些方面的专家......

【问题讨论】:

  • this 有帮助吗?或this?
  • 好的,从这里开始阅读,似乎是 R 中的一些方法来实现这一点。我将阅读更多关于此的内容,但我远未理解这一切。我希望有人可以帮助我,但如果不可能,我可以学习!谢谢!
  • 你可以考虑在gis.stackexchange.com发帖。

标签: r google-maps geolocation latitude-longitude


【解决方案1】:

所以这里发生了几件事。首先,您的数据集似乎具有 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(...) 的滥用,因为重点是多次运行该过程并平均结果,但我只是没有耐心。

【讨论】:

  • 哇...我只是在笑读这个,因为我在第一次阅读时就理解了一半...男人!你是这方面的巫师!我知道只需要采用 depth=0,但我需要将这个“距离”应用于所有数据点......我可以调整它。我可以做的另一件事是在单独的 DF 中提取不同的纬度/经度并在其上运行代码。然后将其用作查找以应用于 2.4mRows...我正在运行一个 8Gig @64bit 的 4 核快速处理器...我希望它会工作。我明天会试试这个并提供反馈。
  • 刚刚数了一下,我有 116k 行不同的经纬度。我将从这个开始。
  • 是的,实际上并行化有很大帮助。查看我的编辑(最后)。
  • 这是一个很好的答案。这是我 2015 年的第一篇笔记。
  • 哇!你真的是个巫师!祝你和家人2015年快乐。我想提一提的是,我从 NOAA 提取的原始数据集在这里引起了很大的讨论:wattsupwiththat.com/2014/12/30/ph-sampling-density我希望加上距离岸边的距离,这将有助于更多的讨论和分析。
猜你喜欢
  • 1970-01-01
  • 2018-06-28
  • 1970-01-01
  • 1970-01-01
  • 2012-07-07
  • 2011-10-01
  • 1970-01-01
  • 1970-01-01
  • 2016-01-11
相关资源
最近更新 更多