2017-02-14 11 views
2

私はオブジェクトに最も近い気象観測所を決定するためにスクリプトをrで書いています。例外なく、私のコードを実行すると、気象ステーション3(テーブルステーションのカムで示される)が自動的に返されます。これがなぜそのようなのでしょうか?私はまた、できるだけインデックスの代わりに駅名を取得したいと思います。オブジェクトに最も近い気象観測所を決定する

コード:

earth.dist <- function (long1, lat1, long2, lat2) 
{ 
rad <- pi/180 
a1 <- lat1 * rad 
a2 <- long1 * rad 
b1 <- lat2 * rad 
b2 <- long2 * rad 
dlon <- b2 - a2 
dlat <- b1 - a1 
a <- (sin(dlat/2))^2 + cos(a1) * cos(b1) * (sin(dlon/2))^2 
c <- 2 * atan2(sqrt(a), sqrt(1 - a)) 
R <- 6378.145 
d <- R * c 
return(d) 
} 


for (i in 1:length(Object$Lat)) 
{ 
    for (j in 1:length(Station$Lat)) 
    { 
     a[j] <- earth.dist(Station$Lon[j], Station$Lat[j], Object$Lon[i], Object$Lat[i]) 
    } 
    index <- which(min(a) %in% a) 
    Object$Station[i] = Station$Station[index] 
} 

駅表:

Station Lat Lon 
SF 37.7749 -122.4194 
CH 41.8781 -87.6298 
Cam 52.2053 -0.1218 

オブジェクト表:我々は、/ R/T地球の距離CALCSワット車輪の再発明をする必要はありません

Object Lat Lon 
1 38.983 -123.092 
2 36.941 -121.767 
3 36.121 -121.084 
4 38.415 -121.787 
5 36.854 -121.362 
6 38.651 -121.218 
7 37.314 -120.386 
8 36.158 -119.8514 
9 38.599 -121.54 
10 35.335 -120.734 
11 34.841 -120.212 
12 38.004 -122.02 
13 37.599 -122.052 
14 38.0903 -122.5267 
15 37.664 -121.885 
16 51.50853 -0.076132 

答えて

1

データ:

read.table(text="Station Lat Lon 
SF 37.7749 -122.4194 
CH 41.8781 -87.6298 
Cam 52.2053 -0.1218", stringsAsFactors=FALSE, header=TRUE) -> stations 

read.table(text="Object Lat Lon 
1 38.983 -123.092 
2 36.941 -121.767 
3 36.121 -121.084 
4 38.415 -121.787 
5 36.854 -121.362 
6 38.651 -121.218 
7 37.314 -120.386 
8 36.158 -119.8514 
9 38.599 -121.54 
10 35.335 -120.734 
11 34.841 -120.212 
12 38.004 -122.02 
13 37.599 -122.052 
14 38.0903 -122.5267 
15 37.664 -121.885 
16 51.50853 -0.076132", stringsAsFactors=FALSE, header=TRUE) -> objs 

コード:

library(geosphere) 
library(tidyverse) 

find_closest_station <- function(lon, lat) { 

    mutate(stations, dist=map2_dbl(Lon, Lat, ~distHaversine(c(lon, lat), c(.x, .y)))) %>% 
    top_n(-1, wt=dist) %>% 
    .$Station 

} 

mutate(objs, wx_st=map2_chr(Lon, Lat, find_closest_station)) 
## Object  Lat   Lon wx_st 
## 1  1 38.98300 -123.092000 SF 
## 2  2 36.94100 -121.767000 SF 
## 3  3 36.12100 -121.084000 SF 
## 4  4 38.41500 -121.787000 SF 
## 5  5 36.85400 -121.362000 SF 
## 6  6 38.65100 -121.218000 SF 
## 7  7 37.31400 -120.386000 SF 
## 8  8 36.15800 -119.851400 SF 
## 9  9 38.59900 -121.540000 SF 
## 10  10 35.33500 -120.734000 SF 
## 11  11 34.84100 -120.212000 SF 
## 12  12 38.00400 -122.020000 SF 
## 13  13 37.59900 -122.052000 SF 
## 14  14 38.09030 -122.526700 SF 
## 15  15 37.66400 -121.885000 SF 
## 16  16 51.50853 -0.076132 Cam 
0

は基本的に、あなたのwhich()ロジックは少しオフになっています。常にTRUEを返すの下に説明するために、そうしかし、以下which()戻り1.

which("R" %in% LETTERS) 
# [1] 1 

TRUE 18日の手紙である:

which(LETTERS == "R") 
# [1] 18 

だから、単純で

index <- which(min(a) %in% a) 

を置き換えます。

index <- which(a == min(a)) 
また、あなたが applyとし、得られインデックス(二つのベクトルを渡す) sapplyでネストされた forを置き換えることができますようにベースRに滞在する必要がある場合は、ソリューションを適用を検討

dist.matrix <- sapply(seq(nrow(Station)), function(x, y) 
          earth.dist(Station$Lon[x], Station$Lat[x], 
             Object$Lon[y], Object$Lat[y]), 
         seq(nrow(Object))) 

Object$Station <- apply(dist.matrix, 1, function(i) Station$Station[which(i == min(i))]) 

Object 
# Object  Lat   Lon  Station 
# 1  1 38.98300 -123.092000   SF 
# 2  2 36.94100 -121.767000   SF 
# 3  3 36.12100 -121.084000   SF 
# 4  4 38.41500 -121.787000   SF 
# 5  5 36.85400 -121.362000   SF 
# 6  6 38.65100 -121.218000   SF 
# 7  7 37.31400 -120.386000   SF 
# 8  8 36.15800 -119.851400   SF 
# 9  9 38.59900 -121.540000   SF 
# 10  10 35.33500 -120.734000   SF 
# 11  11 34.84100 -120.212000   SF 
# 12  12 38.00400 -122.020000   SF 
# 13  13 37.59900 -122.052000   SF 
# 14  14 38.09030 -122.526700   SF 
# 15  15 37.66400 -121.885000   SF 
# 16  16 51.50853 -0.076132   Cam 

でも基地Rのouter列マージンによって後者applyランようsapply上記出力に比べ転置dist.matrixで動作します

dist.mat <- outer(seq(nrow(Station)), seq(nrow(Object)), 
       function(x, y) earth.dist(Station$Lon[x], Station$Lat[x], 
              Object$Lon[y], Object$Lat[y])) 

Object$Station <- apply(dist.mat, 2, function(i) Station$Station[which(i == min(i))]) 
関連する問題