代码之家 › 专栏 › 技术社区 › user2602640

第一次通过时间

user2602640 · 技术社区 · 7 年前

我在用包裹 adehabitatLT 计算首次通过时间。我不太明白它是怎么计算出来的。我假设如果我沿着轨迹计算速度,然后用速度来计算1/速度*半径(乘以2?),我应该得到相同的值。但我没有。fpt是如何计算的,我遗漏了什么?该函数的源代码似乎在C。。。

library(dplyr)
library(adehabitatLT)

radii <- seq(1, 2, 1)

sub <- structure(list(DateTime = structure(c(1497961320, 1497961670, 
1497961833, 1497961975, 1497962075, 1497962211, 1497962417, 1497962584, 
1497962857, 1497963033, 1497963498, 1497963692, 1497963844, 1497963964, 
1497964142, 1497964703, 1497964823, 1497965165, 1497965317, 1497965436, 
1497965610, 1497965776, 1497965923), class = c("POSIXct", "POSIXt"
), tzone = "PST8PDT"), Easting = c(549550.463222071, 549616.446682493, 
549624.93893591, 549621.287381029, 549634.852043099, 549653.472247283, 
549706.355849625, 549752.824737699, 549805.695860071, 549809.742959711, 
549822.329239614, 549850.767478972, 549870.082191359, 549878.454496366, 
549842.209234682, 549810.803726788, 549810.817539391, 549811.560290745, 
549811.995877951, 549785.99096682, 549697.447940334, 549634.423314139, 
549590.138532179), Northing = c(5393845.268479, 5393848.51944054, 
5393886.83599659, 5393933.60504873, 5393964.18478323, 5393997.25495826, 
5394000.39111759, 5393980.0146384, 5393959.69535379, 5394001.19661814, 
5394050.66656145, 5393988.99899494, 5394001.73248796, 5394011.92310602, 
5394057.73545498, 5393997.75985148, 5393996.20363397, 5393995.43205644, 
5393996.10292559, 5393999.20718151, 5393975.85546091, 5393962.73581927, 
5393902.53644517), Burst = c("1", "1", "1", "1", "1", "1", "1", 
"1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", 
"1", "1", "1")), class = c("tbl_df", "tbl", "data.frame"), row.names = c(NA, 
-23L))

计算弹道和fpt

L1.traj <- as.ltraj(data.frame(sub$Easting, sub$Northing), sub$DateTime, 
    id = 1, burst = sub$Burst,
    proj4string = CRS("+proj=utm +zone=11 +ellps=WGS84"))

L1.fpt <- fpt(L1.traj, radii, unit = "seconds")
L1.fpt <- as.data.frame(L1.fpt[[1]])

temp <- ld(L1.traj) %>%
    mutate(DateTime = date) %>%
    dplyr:::select(DateTime, dist, dt) %>%
    # time to pass a radius of 1 m and a diameter of 2 m
    mutate(TimeRadius1 = dt/dist * 1,
        TimeDiameter1 = TimeRadius1 * 2) %>%
    cbind(L1.fpt)

通过半径/直径的时间和fpt的输出之间有一个相当大的差异(而且不一致)。

1 回复 | 直到 7 年前

duckmayr 7 年前

当你需要知道R函数是如何工作的,但它实际上只是调用了用C(或C++或FORTRAN……)编写的函数,你可以通过下载源程序包或者通过查看GITHUB上的源的只读镜像来获得它们的源代码。例如,如果我用谷歌

adehabitatLT cran github

第一个打击是 this “ead是CranR包存储库的唯一镜像。阿得特”。然后你可以在 tests.c fptt 函数调用的 fipati 由 fipatir 从R调用的函数,如果你认为它们的计算有问题,你可以自己判断;他们是这样计算的:

/* compute the FPT for ONE relocation */
void fptt(double *x, double *y, double *t, int pos, double radius, double *fptto, int nlo)
{
    /* Declaration */
    int ok, pos2, naar, naav, na;
    double di, dt, dt2, di2, fptar, fptav;

    ok = 0;
    di = 0;
    di2 = 0;
    dt = 0;
    dt2 = 0;
    naar = 1;
    naav = 1;
    fptar = 0;
    fptav = 0;


    /* Search of the first loc outside the circle (before) */
    pos2 = pos;
    while (ok == 0) {
    pos2 = pos2 - 1;
    if (pos2 > 0) {
        dtmp(x[pos2], x[pos], y[pos2], y[pos], &di);
        if (di >= radius)
        ok = 1;
    } else {
        ok = 1;
        naar = 0;
    }
    }

    /* computes the linear approximation */
    if (naar > 0) {
    dt = fabs(t[pos] - t[pos2]);
    dt2 = fabs(t[pos] - t[(pos2+1)]);
    dtmp(x[(pos2+1)], x[pos], y[(pos2+1)], y[pos], &di2);
    fptar = dt2 + ( (dt - dt2) * (radius - di2) / (di - di2) );
    }


    /* Search of the first loc outside the circle (after) */
    pos2 = pos;
    ok = 0;
    while (ok == 0) {
    pos2 = pos2 + 1;
    if (pos2 <= nlo) {
        dtmp(x[pos2], x[pos], y[pos2], y[pos], &di);
        if (di >= radius)
        ok = 1;
    } else {
        ok = 1;
        naav = 0;
    }
    }

    /* Computes linear approximation */
    if (naav > 0) {
    dt = fabs(t[pos2] - t[pos]);
    dt2 = fabs(t[(pos2-1)] - t[pos]);
    dtmp(x[(pos2-1)], x[pos], y[(pos2-1)], y[pos], &di2);
    fptav = dt2 + ( (dt - dt2) * (radius - di2) / (di - di2) );
    }

    na = naar * naav;
    if (na > 0) {
    *fptto = fptar + fptav;
    } else {
    *fptto = -1;
    }

}