2014-12-29 21 views
6

Tôi đã bắt đầu dự án nguồn mở "miễn phí" để tạo một tập dữ liệu mới cho pH của đại dương trái đất.khoảng cách kinh độ vĩ độ đại dương từ bờ

tôi bắt đầu ngay từ đầu phiên dữ liệu thiết lập từ NOAA và tạo ra một 2,45 triệu hàng dữ liệu thiết lập với các cột:

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

Phương pháp tài liệu HERE.

Tập dữ liệu HERE.

Mục tiêu của tôi bây giờ là "đủ điều kiện" mỗi hàng (2,45m) ... để làm như vậy, tôi cần tính toán khoảng cách từ mỗi điểm Lat/Long đến bờ gần nhất.

Vì vậy, tôi đang tìm kiếm một phương pháp mà sẽ mất Trong: Lat/Long Out: Khoảng cách (km từ bờ)

Với điều này, tôi có thể hội đủ điều kiện nếu điểm dữ liệu có thể bị ảnh hưởng từ ô nhiễm bờ, ví dụ như nước thải của thành phố gần đó.

Tôi đã tìm kiếm một phương pháp để thực hiện việc này, nhưng tất cả dường như cần gói/phần mềm mà tôi không có.

Nếu ai đó sẵn sàng trợ giúp, tôi sẽ đánh giá cao. Hoặc nếu bạn biết một phương pháp dễ dàng (miễn phí) để thực hiện việc này, vui lòng cho tôi biết ...

Tôi có thể làm việc trong lập trình R, Shell script stuff, nhưng không phải chuyên gia ...

+1

[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

+0

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! –

+0

Bạn có thể xem xét việc đăng bài này trên http://gis.stackexchange.com/. – jlhoward

Trả lời

7

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.

+0

Chà ... Tôi chỉ cười khi đọc điều này, bởi vì tôi hiểu một nửa trong số đó là lần đầu đọc ... Đàn ông! Bạn là một thuật sĩ ở đây! Tôi hiểu sự cần thiết phải đi sâu = 0 chỉ, nhưng tôi sẽ cần phải áp dụng "khoảng cách" này cho tất cả các điểm dữ liệu ... Tôi có thể điều chỉnh cho nó. Điều khác tôi có thể làm là trích xuất lat/dài riêng biệt trong một DF riêng biệt và chạy mã trên đó. Sau đó sử dụng nó như là một tra cứu để áp dụng cho 2.4mRows ... Tôi đang chạy một bộ xử lý nhanh 4 lõi với 8Gig @ 64bit ... Tôi hy vọng nó sẽ làm việc. Tôi sẽ cố gắng vào ngày mai và đưa ra phản hồi. –

+0

Chỉ cần đếm, tôi có 116k hàng Lat/Long riêng biệt. Tôi sẽ bắt đầu với điều này. –

+0

Vâng, sự song song thực sự giúp ích rất nhiều. Xem các chỉnh sửa của tôi (ở cuối). – jlhoward

Các vấn đề liên quan