Vì vậy, có một số điều đang xảy ra ở đây. Đầu tiên, tập dữ liệu của bạn dường như có độ pH so với chiều sâu. Vì vậy, trong khi có ~ 2.5MM hàng, chỉ có ~ 200.000 hàng với độ sâu = 0 - vẫn còn rất nhiều.
Thứ hai, để có khoảng cách đến bờ biển gần nhất, bạn cần một shapefile đường bờ biển. May mắn thay điều này có sẵn here, tại số Natural Earth website tuyệt vời. Thứ ba, dữ liệu của bạn dài/lat (vì vậy, đơn vị = độ), nhưng bạn muốn khoảng cách tính bằng km, vì vậy bạn cần phải chuyển đổi dữ liệu của mình (dữ liệu đường bờ biển ở trên cũng dài/lat và cũng cần phải được chuyển đổi). Một vấn đề với biến đổi là dữ liệu của bạn rõ ràng là toàn cầu, và bất kỳ chuyển đổi toàn cục nào cũng sẽ không phải là phẳng. Vì vậy, độ chính xác sẽ phụ thuộc vào vị trí thực tế. Cách đúng đắn để thực hiện điều này là lưới dữ liệu của bạn và sau đó sử dụng một tập hợp các phép biến đổi phẳng phù hợp với bất kỳ điểm lưới nào mà bạn đang ở. Tuy nhiên, điều này nằm ngoài phạm vi của câu hỏi này, vì vậy chúng tôi sẽ sử dụng chuyển đổi toàn cầu (mollweide) chỉ để cung cấp cho bạn một ý tưởng về cách nó được thực hiện trong 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")
Vì vậy, đây đọc dữ liệu của bạn, trích xuất hàng nơi Depth==0
, và chuyển đổi đó để phản đối một SpatialPoints. Sau đó, chúng tôi đọc cơ sở dữ liệu đường bờ biển được tải xuống từ liên kết ở trên thành đối tượng SpatialLines. Sau đó, chúng tôi chuyển đổi cả hai sang phép chiếu Mollweide sử dụng spTransform(...)
, sau đó chúng tôi sử dụng gDistance(...)
trong gói rgeos
để tính khoảng cách tối thiểu giữa mỗi điểm và bờ biển gần nhất.
Một lần nữa, điều quan trọng cần nhớ là mặc dù tất cả các chữ số thập phân, khoảng cách này chỉ xấp xỉ.
Một vấn đề rất lớn là tốc độ: quá trình này mất ~ 2 phút cho 1000 khoảng cách (trên hệ thống của tôi), do đó, để chạy tất cả 200.000 khoảng cách sẽ mất khoảng 6,7 giờ. Một lựa chọn, về mặt lý thuyết, sẽ là tìm một cơ sở dữ liệu bờ biển với độ phân giải thấp hơn.
Mã bên dưới sẽ tính tất cả 201.000 khoảng cách.
## not run
## estimated run time ~ 7 hours
result <- sapply(1:length(sp.points), function(i)gDistance(sp.points[i],coast))
EDIT: Cảm nhận OP về các lõi đã cho tôi để suy nghĩ rằng đây có thể là một trường hợp mà sự cải thiện từ song song có thể xứng đáng với công sức. Vì vậy, đây là cách bạn sẽ chạy này (trên Windows) bằng cách xử lý song song.
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
Sử dụng 4 lõi cải thiện tốc độ xử lý theo yếu tố 3. Vì vậy, kể từ 1000 khoảng cách mất khoảng một phút, 100.000 sẽ mất ít hơn 2 giờ.
Lưu ý rằng việc sử dụng times=1
là lạm dụng microbenchmark(...)
thực sự, vì toàn bộ vấn đề là chạy quá trình nhiều lần và trung bình kết quả, nhưng tôi không có sự kiên nhẫn.
[this] (http://stackoverflow.com/questions/27384403/calculating-minimum-distance-between-a-point-and-the-coast-in-the-uk/27391421#27391421) có giúp được không? hoặc [this] (http://stackoverflow.com/questions/21295302/calculating-minimum-distance-between-a-point-and-the-coast/21302609#21302609)? – jlhoward
Ok đọc từ đây, có vẻ là một số cách trong R để thực hiện điều này. Tôi sẽ đọc thêm về điều này, nhưng tôi không hiểu hết về điều này. Tôi đã hy vọng một người nào đó có thể đưa tay cho tôi, nhưng nếu không thể, tôi có thể học! Cảm ơn! –
Bạn có thể xem xét việc đăng bài này trên http://gis.stackexchange.com/. – jlhoward