fork download 
  1. # 各輸入檔案乃採用原發問者上傳之 csv 和 nc 檔
  2.  
  3. library(magrittr)
  4. library(data.table)
  5. library(ncdf4)
  6.  
  7. # 各種讀資料
  8. analog.indices <-
  9.   fread("HIRAM_WRF_analog.indices(1979-2005).csv",
  10.         col.names = c("day", "i1", "i2")) %>%
  11.   .[1:9862, -1] %>% # 發問者要求刪除 9862 列之後；第 1 欄刪除
  12.   as.matrix # 轉 matrix 之後較好懂
  13. analog.weights <-
  14.   fread("HIRAM_WRF_weights(1979-2005).csv",
  15.         col.names = c("day", "w1", "w2")) %>%
  16.   .[1:9862, -1] %>% # 發問者要求刪除 9862 列之後；第 1 欄刪除
  17.   as.matrix # 轉 matrix 之後較好懂
  18. HIRAM_WRF_data <- nc_open("197901-200512_pr_axis_time_domain.nc")
  19. hiram_wrf_lon <- ncvar_get(HIRAM_WRF_data, "lon")
  20. hiram_wrf_lat <- ncvar_get(HIRAM_WRF_data, "lat")
  21. hiram_wrf_time <- ncvar_get(HIRAM_WRF_data, "time")
  22. hiram_wrf_pr <- ncvar_get(HIRAM_WRF_data, "pr")
  23. dim(hiram_wrf_pr)
  24.  
  25. pr.new <-
  26.   array(
  27.     NA_real_,
  28.     dim = c(dim(hiram_wrf_pr)[1:2], nrow(analog.indices)),
  29.     # 填上合適的 dimnames（尚無實際作用）
  30.     dimnames = list(
  31.       hiram_wrf_lon %>% sprintf("%1.2f", .),
  32.       hiram_wrf_lat %>% sprintf("%1.2f", .),
  33.       1:nrow(analog.indices)
  34.     )
  35.   )
  36.  
  37. ##########################################
  38. # 方法 A：直接列式子
  39. # analog.indices 逐列處理（1:10227）
  40. # 每一圈會為 pr.new 填上一層資料
  41. pr.new.A <- pr.new
  42. for (i in 1:nrow(analog.indices)) {
  43.   pr.new.A[, , i] <-
  44.     hiram_wrf_pr[, , analog.indices[i, 1]] * analog.weights[i, 1] +
  45.     hiram_wrf_pr[, , analog.indices[i, 2]] * analog.weights[i, 2]
  46. }
  47.  
  48. ##########################################
  49. # 方法 B：二層巢式迴圈
  50. # 外層 i 照舊走直的，內層 j 走橫的
  51. # 每一小圈會計算部份累積的總和
  52. # 每一大圈會為 pr.new 填上一層資料
  53. pr.new.B <- pr.new
  54. for (i in 1:nrow(analog.indices)) {
  55.   partialSum <- 0.0
  56.   for (j in 1:ncol(analog.indices)) {
  57.     partialSum <-
  58.       partialSum +
  59.       hiram_wrf_pr[, , analog.indices[i, j]] * analog.weights[i, j]
  60.   }
  61.   pr.new.B[, , i] <- partialSum
  62.   # 以上用了一個以迴圈求累積加總的常見小技巧
  63.   # 沒看過的話，看看這個求 4+6+7 的例子：
  64.   # a <- c(4,6,7); b <- 0; for(i in 1:3){b <- b + a[i]}; print(b)
  65. }
  66.  
  67. ###################################################
  68. # 檢查第11日第 (1,1) 格的結果
  69. analog.indices[11,]
  70. # i1 i2 
  71. # 10 11 
  72. analog.weights[11,]
  73. #        w1        w2 
  74. # 0.2276517 0.7622840 
  75. print(hiram_wrf_pr[1, 1, 10] * 0.2276517 + hiram_wrf_pr[1, 1, 11] * 0.762284) # 0.007988894
  76. print(pr.new.A[1, 1, 11]) # 0.007988894
  77. print(pr.new.B[1, 1, 11]) # 0.007988894
  78. identical(pr.new.A, pr.new.B) # 二者全等
  79.  
Success #stdin #stdout #stderr 0.26s 188800KB
comments (?)
stdin 
Standard input is empty
stdout 
Standard output is empty
stderr 
Error in library(ncdf4) : there is no package called ‘ncdf4’
Execution halted
https://ideone.com/nRizTB
language:
R (R 3.5.2)
created:
5 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!