r - 确定最近的站点并从该位置选择另一个变量

Question

我有收集数据的研究地点和附近的气象站，其中包含有关温度和降水的信息。 我想将我在研究地点的每日数据与最近的气象站的天气信息配对。 我认为，要做到这一点，我需要一个两步过程，首先选择离研究地点最近的气象站，然后用天气数据创建一个新变量。

这是我的数据快照：

# study sites
site <- rep(LETTERS[1:3], 5)
siteLat <- rep(c(41, 42, 44), 5)
siteLon <- rep(c(68, 62, 63), 5)
siteDate <- rep(1:5, 3)
dfSites <- data.frame(cbind(site, siteLat, siteLon, siteDate))
# weather stations
station <- rep(letters[1:3], 5)
stationLat <- rep(c(40, 43, 45), 5)
stationLon <- rep(c(67, 61, 64), 5)
stationDate <- rep(1:5, 3)
temp <- sample(10:20, 15, replace=TRUE)
dfStation <- data.frame(cbind(station, stationLat, stationLon, stationDate, temp))

我试图用这条线来确定哪个车站最近，但我只得到一排距离。

distVincentyEllipsoid(df2[c("recvLon", "recvLat")], weather[c("lon", "lat")])

计算完所有距离后，我有点不确定下一步，但我认为我需要一些东西来选择最近的车站和匹配日期。这是我想出的最好的：

dfSites %>% 
mutate(closestStation = ???,
temp1 = temp[station == closestStation & stationDate == siteDate])

最终结果是我的研究站点数据框，其中包含来自最近气象站的附加温度列。

Answer 1

我认为distVincentyEllipsoid(p1, p2, ...)试图找到p1的第一点与p2的第一点、p1的第二点与p2的第二点之间的距离，等等。你需要的是按照 *"p1 中第一个针对所有p2，第二个在p1中针对所有p2，等等)进行扩展。

调整您的代码以调用dfSites和dfStation(而不是df2/weather)，以下内容应该适合您。(我将删除其中一个带有dfStation[-1,...]的站点，只是为了清楚地识别哪个维度代表站点与站点。

alldists <- sapply(seq_len(nrow(dfSites)), function(i) {
distVincentyEllipsoid(dfSites[i,c("siteLon","siteLat")],
dfStation[-1,c("stationLon","stationLat")])
})
alldists
#           [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]
#  [1,] 786180.9 123505.1 228960.0 786180.9 123505.1 228960.0 786180.9 123505.1
#  [2,] 481351.6 269760.4 122086.2 481351.6 269760.4 122086.2 481351.6 269760.4
#  [3,] 119427.7 565573.7 484015.5 119427.7 565573.7 484015.5 119427.7 565573.7
#  [4,] 786180.9 123505.1 228960.0 786180.9 123505.1 228960.0 786180.9 123505.1
#  [5,] 481351.6 269760.4 122086.2 481351.6 269760.4 122086.2 481351.6 269760.4
#  [6,] 119427.7 565573.7 484015.5 119427.7 565573.7 484015.5 119427.7 565573.7
#  [7,] 786180.9 123505.1 228960.0 786180.9 123505.1 228960.0 786180.9 123505.1
#  [8,] 481351.6 269760.4 122086.2 481351.6 269760.4 122086.2 481351.6 269760.4
#  [9,] 119427.7 565573.7 484015.5 119427.7 565573.7 484015.5 119427.7 565573.7
# [10,] 786180.9 123505.1 228960.0 786180.9 123505.1 228960.0 786180.9 123505.1
# [11,] 481351.6 269760.4 122086.2 481351.6 269760.4 122086.2 481351.6 269760.4
# [12,] 119427.7 565573.7 484015.5 119427.7 565573.7 484015.5 119427.7 565573.7
# [13,] 786180.9 123505.1 228960.0 786180.9 123505.1 228960.0 786180.9 123505.1
# [14,] 481351.6 269760.4 122086.2 481351.6 269760.4 122086.2 481351.6 269760.4
#           [,9]    [,10]    [,11]    [,12]    [,13]    [,14]    [,15]
#  [1,] 228960.0 786180.9 123505.1 228960.0 786180.9 123505.1 228960.0
#  [2,] 122086.2 481351.6 269760.4 122086.2 481351.6 269760.4 122086.2
#  [3,] 484015.5 119427.7 565573.7 484015.5 119427.7 565573.7 484015.5
#  [4,] 228960.0 786180.9 123505.1 228960.0 786180.9 123505.1 228960.0
#  [5,] 122086.2 481351.6 269760.4 122086.2 481351.6 269760.4 122086.2
#  [6,] 484015.5 119427.7 565573.7 484015.5 119427.7 565573.7 484015.5
#  [7,] 228960.0 786180.9 123505.1 228960.0 786180.9 123505.1 228960.0
#  [8,] 122086.2 481351.6 269760.4 122086.2 481351.6 269760.4 122086.2
#  [9,] 484015.5 119427.7 565573.7 484015.5 119427.7 565573.7 484015.5
# [10,] 228960.0 786180.9 123505.1 228960.0 786180.9 123505.1 228960.0
# [11,] 122086.2 481351.6 269760.4 122086.2 481351.6 269760.4 122086.2
# [12,] 484015.5 119427.7 565573.7 484015.5 119427.7 565573.7 484015.5
# [13,] 228960.0 786180.9 123505.1 228960.0 786180.9 123505.1 228960.0
# [14,] 122086.2 481351.6 269760.4 122086.2 481351.6 269760.4 122086.2

(因为我们有 14 行，所以每行都是您的一个站点。您不应该执行[-1,]索引，只需知道哪一行/列即可。由此我们知道，站点A和站点b之间的差异为481351.6米(第一列，第二排)。

从这里，只需找到列最小值：

apply(alldists, 2, which.min)
#  [1] 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2

建议离站点A最近的车站是b(which.min将返回第一个最小值，它不表示连接)。

现在，dfStation[apply(alldists, 2, which.min),]为您提供了 15 行台站数据，可以轻松cbind或以其他方式与dfSites相结合。

---

dplyr选项：

dfSites %>%
mutate(
station_i = purrr::map2_int(
siteLat, siteLon,
~ which.min(geosphere::distVincentyEllipsoid(
cbind(.x,.y), dfStation[-1,c("stationLon","stationLat")]))
),
station = as.character(dfStation$station)[ station_i ]
)
#    site siteLat siteLon siteDate station_i station
# 1     A      41      68        1         3       c
# 2     B      42      62        2         1       a
# 3     C      44      63        3         2       b
# 4     A      41      68        4         3       c
# 5     B      42      62        5         1       a
# 6     C      44      63        1         2       b
# 7     A      41      68        2         3       c
# 8     B      42      62        3         1       a
# 9     C      44      63        4         2       b
# 10    A      41      68        5         3       c
# 11    B      42      62        1         1       a
# 12    C      44      63        2         2       b
# 13    A      41      68        3         3       c
# 14    B      42      62        4         1       a
# 15    C      44      63        5         2       b

---

通过它们的外表产品可以看到轻微(10-15%)的速度提高。

outer(seq_len(nrow(dfSites)), seq_len(nrow(dfStation)),
function(i,j) geosphere::distVincentyEllipsoid(dfSites[i,2:3], dfStation[j,2:3]))

这还会返回一个mxn矩阵(工作站行)，然后您apply(...)该矩阵以获取最接近的索引。(我希望有更大的性能提升，因为distVincentyEllipsoid只调用一次......