fix rdist.earth

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
+}

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
 }

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
-}

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

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

R/tools.R

@@ -59,3 +59,11 @@ mode <- function(x) {
     names() %>%
     as.numeric()
 }
+
+`%||%` <- function(x, y) {
+  if (is.null(x)) {
+    y
+  } else {
+    x
+  }
+}

scripts/met840/Step03_缺失信息检测与插值.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()

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)
 })