diff --git a/R/find_near.R b/R/find_near.R new file mode 100644 index 0000000..7fe9648 --- /dev/null +++ b/R/find_near.R @@ -0,0 +1,13 @@ +#' find_near +#' @param p a vector of length 2, `[lon, lat]` +#' @export +find_near <- function(p, st) { + loc <- st[, .(lon, lat)] + loc %<>% + check_matrix() + dist <- rdist.earth2(loc, p) + i <- which.min(dist) + s <- st[i, ] |> + cbind(dist = dist[i]) + s +} diff --git a/R/interp_approx.R b/R/interp_approx.R index e1cf9d3..200a327 100644 --- a/R/interp_approx.R +++ b/R/interp_approx.R @@ -21,19 +21,18 @@ na_approx.dtime <- function(x, maxgap = 5, ...) { #' @param .parallel boolean #' #' @export -interp_approx <- function( - xx, stationInfo, maxgap = 5, verbose = TRUE, - .parallel = FALSE, ...) { +interp_approx <- function(xx, stationInfo, maxgap = 5, verbose = TRUE, + .parallel = FALSE, ...) +{ I <- with(stationInfo, which(n_miss > 0 & gap_min <= maxgap)) # gap satisfied N <- length(I) step <- floor(N / 10) - FUN <- ifelse(.parallel, `%dopar%`, `%do%`) - xx[I] <- FUN(foreach(x = xx[I], i = icount()), { - if (verbose) { - runningId(i, step, N, "interp_approx") - } + `%dof%` <- ifelse(.parallel, `%dopar%`, `%do%`) + + xx[I] <- foreach(x = xx[I], i = icount()) %dof% { + verbose && runningId(i, step, N, "interp_approx") na_approx.dtime(x, maxgap) - }) - return(xx) + } + xx } diff --git a/R/rdist.earth.R b/R/rdist.earth.R index c2e0ee7..0320e4c 100644 --- a/R/rdist.earth.R +++ b/R/rdist.earth.R @@ -10,38 +10,13 @@ #' } #' #' @export -rdist.earth <- function(x1, x2 = x1) { - R <- 6378.388 - # coslat1 <- cos(x1[, 2]) sinlat1 <- sin(x1[, 2]) coslon1 <- cos(x1[, 1]) sinlon1 <- sin(x1[, 1]) - coslat1 <- cos((x1[, 1] * pi) / 180) - sinlat1 <- sin((x1[, 1] * pi) / 180) - coslon1 <- cos((x1[, 2] * pi) / 180) - sinlon1 <- sin((x1[, 2] * pi) / 180) - - coslat2 <- cos((x2[, 1] * pi) / 180) - sinlat2 <- sin((x2[, 1] * pi) / 180) - coslon2 <- cos((x2[, 2] * pi) / 180) - sinlon2 <- sin((x2[, 2] * pi) / 180) - pp <- cbind(coslat1 * coslon1, coslat1 * sinlon1, sinlat1) %*% t(cbind( - coslat2 * coslon2, coslat2 * sinlon2, - sinlat2 - )) - return(R * acos(ifelse(abs(pp) > 1, 1 * sign(pp), pp))) -} +rdist.earth <- function(x1, x2 = NULL) { + if (is.vector(x1)) x1 %<>% t() + if (is.vector(x2)) x2 %<>% t() + x2 = x2 %||% x1 -rdist.earth2 <- function(x1, x2 = NULL) { - if (is.vector(x1)) { - x1 %<>% - t() - } - if (is.vector(x2)) { - x2 %<>% - t() - } - x1 %<>% - check_matrix() - x2 %<>% - check_matrix() + x1 %<>% check_matrix() + x2 %<>% check_matrix() # miles = FALSE, R = NULL R = R %||% ifelse(miles, 3963.34, 6378.388) R <- 6378.388 @@ -76,17 +51,3 @@ check_matrix <- function(x) { } x } - -#' find_near -#' @param p a vector of length 2, `[lon, lat]` -#' @export -find_near <- function(p, st) { - loc <- st[, .(lon, lat)] - loc %<>% - check_matrix() - dist <- rdist.earth2(loc, p) - i <- which.min(dist) - s <- st[i, ] |> - cbind(dist = dist[i]) - s -} diff --git a/R/st-get_moveInfo.R b/R/st-get_moveInfo.R index 8a53fb6..805e0b5 100644 --- a/R/st-get_moveInfo.R +++ b/R/st-get_moveInfo.R @@ -88,7 +88,7 @@ distToCentralPeriod <- function(lon, lat, n_period) { i <- which.max(n_period) # if (nrow(P) > 1) { - dist <- rdist.earth(P[i, , drop = FALSE], P[, , drop = FALSE])[1, ] + dist <- rdist.earth(P[i, , drop = FALSE], P[, , drop = FALSE]) round(dist, 2) } else { 0 diff --git a/R/st-revise_locaiton_multi.R b/R/st-revise_locaiton_multi.R index f8c9ff8..48595be 100644 --- a/R/st-revise_locaiton_multi.R +++ b/R/st-revise_locaiton_multi.R @@ -7,7 +7,7 @@ distToCentralPeriod <- function(lon, lat, n_period) { i <- which.max(n_period) # if (nrow(P) > 1) { - dist <- rdist.earth(P[i, , drop = FALSE], P[, , drop = FALSE])[1, ] + dist <- rdist.earth(P[i, , drop = FALSE], P[, , drop = FALSE]) round(dist, 2) } else { 0 diff --git a/R/tools.R b/R/tools.R index d7642bb..350e851 100644 --- a/R/tools.R +++ b/R/tools.R @@ -59,3 +59,11 @@ mode <- function(x) { names() %>% as.numeric() } + +`%||%` <- function(x, y) { + if (is.null(x)) { + y + } else { + x + } +} diff --git a/man/find_near.Rd b/man/find_near.Rd index ff833d8..df3fc35 100644 --- a/man/find_near.Rd +++ b/man/find_near.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/rdist.earth.R +% Please edit documentation in R/find_near.R \name{find_near} \alias{find_near} \title{find_near} diff --git a/man/rdist.earth.Rd b/man/rdist.earth.Rd index 4ac23de..d2c2bf0 100644 --- a/man/rdist.earth.Rd +++ b/man/rdist.earth.Rd @@ -4,7 +4,7 @@ \alias{rdist.earth} \title{rdist.earth} \usage{ -rdist.earth(x1, x2 = x1) +rdist.earth(x1, x2 = NULL) } \arguments{ \item{x1}{2 columns matrix, [lat, lon]} diff --git "a/scripts/met840/Step03_\347\274\272\345\244\261\344\277\241\346\201\257\346\243\200\346\265\213\344\270\216\346\217\222\345\200\274.R" "b/scripts/met840/Step03_\347\274\272\345\244\261\344\277\241\346\201\257\346\243\200\346\265\213\344\270\216\346\217\222\345\200\274.R" index 372abf6..c31bf6c 100644 --- "a/scripts/met840/Step03_\347\274\272\345\244\261\344\277\241\346\201\257\346\243\200\346\265\213\344\270\216\346\217\222\345\200\274.R" +++ "b/scripts/met840/Step03_\347\274\272\345\244\261\344\277\241\346\201\257\346\243\200\346\265\213\344\270\216\346\217\222\345\200\274.R" @@ -43,7 +43,7 @@ if (file.exists(infile)){ fstation <- "说明文档/中国气候日志数据V3station839.xlsx" stations <- read.xlsx(fstation) -dist <- rdist.earth(stations[, 3:2], miles = F)#lon, lat +dist <- rdist.earth(stations[, 3:2])#lon, lat ## ---------------------- prepare for parallel compute ------------------- site <- stations$StationNo time_day <- as.Date("1951-01-01"):as.Date("2016-06-30") %>% as.Date() diff --git a/tests/testthat/test-interp_main.R b/tests/testthat/test-interp_main.R index aae22b1..8a88dd2 100644 --- a/tests/testthat/test-interp_main.R +++ b/tests/testthat/test-interp_main.R @@ -1,10 +1,10 @@ test_that("interp_main and dtime2mat works", { - data("prcp") - data("st840") + data("prcp") + data("st840") - r <- interp_main(prcp, st840, smax = 200, verbose = FALSE) + r <- interp_main(prcp, st840, smax = 200, verbose = FALSE) - df2 <- dtime2mat(r$x) - expect_equal(ncol(df2), ncol(prcp)) - expect_false(all.equal(df2, prcp) == TRUE) + df2 <- dtime2mat(r$x) + expect_equal(ncol(df2), ncol(prcp)) + expect_false(all.equal(df2, prcp) == TRUE) })