# 各輸入檔案乃採用原發問者上傳之 csv 和 nc 檔

library(magrittr)
library(data.table)
library(ncdf4)

# 各種讀資料
analog.indices <-
  fread("HIRAM_WRF_analog.indices(1979-2005).csv",
        col.names = c("day", "i1", "i2"))
analog.weights <-
  fread("HIRAM_WRF_weights(1979-2005).csv",
        col.names = c("day", "w1", "w2"))
HIRAM_WRF_data <- nc_open("197901-200512_pr_axis_time_domain.nc")
hiram_wrf_lon <- ncvar_get(HIRAM_WRF_data, "lon")
hiram_wrf_lat <- ncvar_get(HIRAM_WRF_data, "lat")
hiram_wrf_time <- ncvar_get(HIRAM_WRF_data, "time")
hiram_wrf_pr <- ncvar_get(HIRAM_WRF_data, "pr")

# 創新的三維 array
pr.new <-
  array(
    NA_real_,
    dim = c(dim(hiram_wrf_pr)[1:2], nrow(analog.indices)),
    # 填上合適的 dimnames（尚無實際作用）
    dimnames = list(
      hiram_wrf_lon %>% sprintf("%1.2f", .),
      hiram_wrf_lat %>% sprintf("%1.2f", .),
      1:nrow(analog.indices)
    )
  )
dim(pr.new) #[1]    41    77 10227

# analog.indices 逐列處理（1:10227）
# 每一圈會為 pr.new 填上一層資料，共 10227 層
for (i in 1:nrow(analog.indices)) {
  pr.new[, , i] <-
    hiram_wrf_pr[, , analog.indices$i1[i]] * analog.weights$w1[i] +
    hiram_wrf_pr[, , analog.indices$i2[i]] * analog.weights$w2[i]
}

# 檢查 pr.new[1,1,11]
analog.indices[11,]
#    day i1 i2
# 1:  10 11 11
analog.weights[11,]
#    day        w1       w2
# 1:  11 0.2276517 0.762284
print(hiram_wrf_pr[1, 1, 10] * 0.2276517 + hiram_wrf_pr[1, 1, 11] * 0.762284) # 0.007988894
print(pr.new[1, 1, 11]) # 0.007988894