根据船速和上次捕获的时间判断鱼的位置是否错误

Question

我有一系列在不同日期时间和不同行程中从船上捕获的鱼的坐标。我如何根据同一行程中最后一条鱼被捕获的时间和假定的船速（比如 10 公里/小时）来确定一条鱼的坐标是否可能不正确（例如，由于转录错误）。

这是一个简单的示例数据集，其中包含 2 次旅行和每次旅行两条鱼。

library(sf)
library(ggplot2)
library(dplyr)
library(lubridate)

datetime <- ymd_hms('2017-05-13 14:00:00', tz = "Etc/GMT+8")
df <- data_frame(DateTimeCapture = c(datetime, datetime + minutes(35), datetime + days(2), 
                                     datetime + days(2) + minutes(20)),
                 Trip = c('1', '1', '2', '2'),
                 Order = c(1, 2, 1, 2),
                 X = c(648635, 648700, 647778, 658889),
                 Y = c(5853151, 5853200, 5854292, 5870000))

# if you prefer to work in sf
df_sf <-  st_as_sf(df, coords = c('X', 'Y'), crs = 32610)

# quick plot
ggplot() + 
  geom_point(data = df, aes(x = X, y = Y, color = Trip))

第二趟两条鱼的距离是19km：

st_distance(df_sf[3:4, ])
Units: m
         [,1]     [,2]
[1,]     0.00 19240.47
[2,] 19240.47     0.00

一艘船不可能在20分钟内行驶19公里。因此，这应该被标记为可能的错误。

我更喜欢使用 sf 的解决方案，但也可以接受使用 sp 的解决方案。它必须是基于 r 的解决方案。

Answer 1

这可能会解决您的问题：

fun1 <- function(k){
  dat <- st_as_sf(df[which(df$Trip == k),], coords = c('X', 'Y'), crs = 32610)
  times <- as.numeric(diff(dat$DateTimeCapture))
  distances <- st_distance(dat)
  distances <- diag(distances[-1,])

  tresh <- 10000/60 # 10km/h is our treshold here
  problematic <- as.numeric(distances/times) > tresh

  if(length(which(problematic)) >= 1){
    v <- matrix(F, nrow = length(dat$Trip))
    v[which(problematic)+1] <- T
    return(v)
  }
  if(length(which(problematic)) == 0){
    v <- matrix(F, nrow = length(dat$Trip))
    return(v)
  }
} # brief explanations below

我的输出

unlist(sapply(unique(df$Trip), fun1, simplify = F))
   11    12    21    22 
FALSE FALSE FALSE  TRUE 

# and now cbinding it into the data frame:
> newcol <- unlist(sapply(unique(df$Trip), fun1, simplify = F))
> df <- cbind(df, newcol)
> df
       DateTimeCapture Trip Order      X       Y newcol
11 2017-05-14 00:00:00    1     1 648635 5853151  FALSE
12 2017-05-14 00:35:00    1     2 648700 5853200  FALSE
21 2017-05-16 00:00:00    2     1 647778 5854292  FALSE
22 2017-05-16 00:20:00    2     2 658889 5870000   TRUE

简要说明

以上 function 检查给定行程是否包含异常。

它计算时间差（times）和距离矩阵（distances)。
现在考虑distances的子或super-diagonal就够了。实际上，对于给定的行程，这些对角线都包含两次连续捕获之间的所有距离。
那么，剩下要做的就是检查是否distance/time > tresh（这里是10km/h）。

现在，function 可以进行调整、润色等。例如您可能希望将 tresh 作为参数传递给函数并使用 missing() 为其指定默认值。

免责声明我稍微编辑了你的数据（在行程 2 中添加了第三个点以获得更有趣的测试用例）：

df <- data.frame(DateTimeCapture = c(datetime, datetime + minutes(35), datetime + days(2), 
                                 datetime + days(2) + minutes(20), datetime + days(2) + minutes(45)),
             Trip = c('1', '1', '2', '2', '2'),
             Order = c(1, 2, 1, 2, 3),
             X = c(648635, 648700, 647778, 658889, 658999),
             Y = c(5853151, 5853200, 5854292, 5870000, 5890978))

Answer 2

sf::st_distance() 生成所有几何之间的距离矩阵。

从这个矩阵中我们可以只提取我们关心的距离，然后使用这些距离来计算行进的速度，如果超过某个阈值则添加一个flag

library(dplyr)

max_speed <- 10 ## km/h


df_sf %>%
    mutate(distance = {
        dist_mat <- sf::st_distance(.)
        distances <- dist_mat[ upper.tri(dist_mat) ]
        idx <- cumsum(2:ncol(dist_mat) - 1)
        distances <- c(0, distances[ idx ] )
        distances[.$Order == 1] <- 0         ## first trip gets 0 distance
        distances
    }) %>%
    mutate( time = as.numeric(difftime(DateTimeCapture, lag(DateTimeCapture))),
                    speed = distance / time) %>%
    mutate( error_flag = speed > max_speed ) 


# 
# Simple feature collection with 4 features and 7 fields
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: 647778 ymin: 5853151 xmax: 658889 ymax: 5870000
# epsg (SRID):    32610
# proj4string:    +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
# # A tibble: 4 x 8
#    DateTimeCapture     Trip  Order distance   time   speed error_flag               geometry
#    <dttm>              <chr> <dbl>    <dbl>  <dbl>   <dbl> <lgl>           <sf_geometry [m]>
# 1  2017-05-14 08:00:00 1      1.00      0     NA    NA     NA         POINT (648635 5853151)
# 2  2017-05-14 08:35:00 1      2.00     81.4   35.0   2.33  F          POINT (648700 5853200)
# 3  2017-05-16 08:00:00 2      1.00      0   2845     0     F          POINT (647778 5854292)
# 4  2017-05-16 08:20:00 2      2.00  19240     20.0 962     T          POINT (658889 5870000)

详情

有关获取距离的第一个 mutate 调用中发生的事情的一些详细信息。

st_distance()函数给出了每个几何体之间的距离矩阵。

dist_mat <- sf::st_distance(df_sf)
dist_mat
# Units: m
#             [,1]        [,2]      [,3]     [,4]
# [1,]     0.00000    81.40025  1427.000 19723.93
# [2,]    81.40025     0.00000  1429.177 19648.30
# [3,]  1427.00035  1429.17739     0.000 19240.47
# [4,] 19723.92752 19648.30072 19240.467     0.00

从这个矩阵中，我们需要 [1, 2]、[2, 3] 和 [3, 4]

处的值

所以首先我们可以采取 upper-triangle

distances <- dist_mat[ upper.tri(dist_mat) ]
distances
# Units: m
# [1]    81.40025  1427.00035  1429.17739 19723.92752 19648.30072 19240.46738

然后取这个向量的第1、3、6个索引

idx <- c(cumsum(2:ncol(dist_mat) - 1))
idx
# [1] 1 3 6

给我们距离

c(0, distances[ idx ] )
# [1]     0.00000    81.40025  1429.17739 19240.46738

根据船速和上次捕获的时间判断鱼的位置是否错误

determine whether location of fish is wrong based on boat speed and time of last capture

gis

r

sp

sf

详情