Re: [問題] 月平均資料

[andrew43 (討厭有好心推文後刪文者)]

# 已先自 http://0rz.tw/JI056 下載 trmm_2010.nc

library(ncdf4)
library(data.table)

obs <- nc_open("trmm_2010.nc")
lon <- ncvar_get(obs, "lon") # 經度切1440份
lat <- ncvar_get(obs, "lat") # 緯度切400份
time <- ncvar_get(obs, "time") # 單位為小時，共365個
precip <- ncvar_get(obs, "r") # 第一維為lon，第二維為lat，笫三維為time

time2 <-
  as.Date(time / 24, format = "%Y-%m-%d", origin = "2010-01-01")
# 先把單位由小時換算成天

##############

# 以每日求12個月平均
# precip.ave.monthly[1, , ] 到 precip.ave.monthly[12, , ] 表示各月均雨量
# 這步很慢，有需要再改寫
precip.ave.monthly <-
  apply(precip, c(1, 2), function(x) {
    tapply(x, month(time2), mean)
  })
image(lon, lat, precip.ave.monthly[6, ,]) # 六月份

# 以每日求全年總平均
precip.ave.overall <- apply(precip, c(1, 2), mean)
image(lon, lat, precip.ave.overall)


par(mfrow = c(3, 4),
    mar = c(2, 2, 2, 2),
    cex = 8 / 12)
for (i in 1:12) {
  image(
    lon,
    lat,
    precip.ave.monthly[i, ,],
    col = gray(0:255 / 256),
    xlab = "",
    ylab = ""
  )
  title(paste("month:", i))
}


※ 引述《AndrewShi (沒有妳的我)》之銘言：
: [問題類型]:
: 程式諮詢(我想用R 做某件事情，但是我不知道要怎麼用R 寫出來)
: [軟體熟悉度]:
: 入門(寫過其他程式，只是對語法不熟悉)
: [問題敘述]:
: 各位大大好，
: 我放入的這筆資料是2010年全球每天的降雨(量)資料，現在我想把每日的降雨量計算成月
: 平均.年平均降雨量，下面我所想到的迴圈是可以畫得出圖來，但畫出來感覺不太正確，所以想請
: 教大大們我的迴圈是否有問題，能否給我一些提點，謝謝。
: p.s:原本的資料型態中降雨值的維度只包含經度和緯度(2維)，所以我用rbind把時間的維
: 度也併到降雨值裡。
: [程式範例]:
: rm(list=ls())
: library(ncdf4)
: TRMM_data <- "C:\\Users\\TOM\\Desktop\\R(資料庫)\\TRMM資料\\trmm_2010.nc"
: obs <- nc_open(TRMM_data)
: print(obs)
: lon <- ncvar_get(obs, "lon")
: lat <- ncvar_get(obs, "lat")
: time <- ncvar_get(obs, "time")
: precip <- ncvar_get(obs,"r")
: time <- matrix(seq(as.Date("2010-01-01"), as.Date("2010-12-31"),1))
: rbind(dim(time),precip[[3]])
: time <- c()
: for(time in seq_along(1:31)){
:   mean(precip)
: }
: time <- c()
: for(time in seq_along(1:365)){
:   mean(precip)
: }
: lon <- lon-180
: #lat <- rev(lat)
: precip <- precip[,,time]
: library(RColorBrewer)
: image(lon,lat,precip,col=rev(brewer.pal(10,"RdBu")))
: library(maptools)
: gpclibPermit()
: data(wrld_simpl)
: plot(wrld_simpl,add=TRUE)
: [環境敘述]:
: [關鍵字]:
: 月平均   nc檔   降雨

--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 60.248.222.1
※ 文章網址: https://www.ptt.cc/bbs/R_Language/M.1540258357.A.42A.html
推 AndrewShi: 非常感謝andrew大，願意花這麼多時間幫我解答，我今天
→ AndrewShi: 也會好好研究程式碼，如果對於你寫的程式碼有疑惑的話
→ AndrewShi: 方便能再請教你嗎XD??
推 AndrewShi: abdrew大~我想請教apply裡的c(1,2)是指把每列.每行(每
→ AndrewShi: 個經度.緯度)的降雨值都帶到function裡面的意思對嗎?!
precip 是一個三維陣列，維度按順序分別是經度、緯度和時間。
apply(precip, c(1, 2), mean) 的意思是，
對 precip 做平均，方式是同經度同緯度而不同時間的值取成一組資料的平均。
這裡的 c(1, 2) 就是指第一和第二個維度，即同經度同緯度的意思。
因此，其結果就由365「層」平均成 1 層。

→ AndrewShi: 另外想請問你是怎麼把時間(月)放到precip.ave.monthly
→ AndrewShi: 的第一個維度裡的呢??
上面如果可以理解，那麼
precip.ave.monthly <-
  apply(precip, c(1, 2), function(x) {
    tapply(x, month(time2), mean)
  })
這句也是類似的。
不過，這次希望同經度同緯度的值在平均時要再按 month(time2) 為組別分別平均。
所以，每個地點會取出 12 個平均，最後結果就是有 12 層。

我這裡是利用 month(time2) 將年月日取到月份。
如果你的資料不只一年，那 month(time2) 要改成 substr(time2, 1, 7) 才行，
不然不能區別同月不同年的情況。
※ 編輯: andrew43 (60.248.222.1), 10/23/2018 17:43:50
※ 編輯: andrew43 (60.248.222.1), 10/23/2018 18:31:26
推 AndrewShi: 感謝andrew大詳細的解答，上面的敘述我能理解，不過我
→ AndrewShi: 很好奇為什麼在單用apply算降雨年平均的時候時間是在
→ AndrewShi: 降雨的第三個維度，而在算月平均的時候，時間變成降雨(
→ AndrewShi: precip.ave.monthly)的第一個維度，而經.緯度則變成第
→ AndrewShi: 二.三個維度，我看不出來哪一段程式碼是在做這樣的處理
→ AndrewShi: 。
→ andrew43: 這應該是因為tapply()一次有12個值，所以apply自動把它
→ andrew43: 安排在第一個維度。沒有哪段碼刻意要求達到這種結果。
→ andrew43: 再接一個aperm(..., c(2,3,1))轉換過來就好了
→ andrew43: 另外，單用apply算降雨年平均後只有二維度，並沒有時間.
推 AndrewShi: 了解，真的非常感謝andrew大，讓我學到很多~
推 AndrewShi: andrew大~想再請教你如果想用迴圈的概念來寫的話，我想
→ AndrewShi: 到的寫法是:for(i in (1:1440)){
→ AndrewShi: for(j in (1:400)){
→ AndrewShi: mean(precip[i,j,time=(1:31)])}}
→ AndrewShi: 但跑出來也只有一個值，所以想請教你我的迴圈寫法(概
→ AndrewShi: 念)是哪裡有出錯嗎??
→ andrew43: 這樣是各別算出每個地點一月份的平均雨量，一次算一個。
→ andrew43: 那還不如寫成apply(precip[,,1:31], c(1,2), mean)
→ andrew43: 你若想用for loop一次算一個地點，要把結果填到一個容器
→ andrew43: 中，例如一個400*1440矩陣中。
推 AndrewShi: 了解，非常感謝~
推 AndrewShi: andrew大~不好意思，想再請教你一個基本問題，把結果
→ AndrewShi: 填到一個矩陣裡，matrix(mean(precip[i,j,time=(1:31),
→ AndrewShi: 1440,400))這樣寫對嗎??
不對。
你似乎不太熟悉一些矩陣或陣列的操作。
請研究以下例子。

a <- array(1:24, c(2, 3, 4))
m <- matrix(NA_real_, nrow = dim(a)[1], ncol = dim(a)[2])
for (i in 1:dim(a)[1]) {
  for (j in 1:dim(a)[2]) {
    m[i, j] <- mean(a[i, j, ])
  }
}

print(a)
print(m)
apply(a, c(1, 2), mean)

※ 編輯: andrew43 (60.248.222.1), 10/26/2018 18:25:27
※ 編輯: andrew43 (60.248.222.1), 10/26/2018 18:27:53
推 AndrewShi: 我懂了，首先先創一個數字由1~24，2列3行共4個(矩陣)
→ AndrewShi: 的陣列a，之後再創一個矩陣m，列和行的數目和陣列a的
→ AndrewShi: 第一和第二個維度一樣(2列3行)，最後再把陣列a每個列和
→ AndrewShi: 行各別的值(共4個)相加取平均後放到矩陣m裡(如:(1+7+13
→ AndrewShi: +19)/4=10)，我這樣的理解應該沒有錯吧 :)?!
推 AndrewShi: 而我的程式應該改為:precip1 <- matrix(precip,1440,40
→ AndrewShi: 0)
→ AndrewShi: for(i in (1:1440)){
→ AndrewShi: for(j in 1:400)){
→ AndrewShi: precip1[i,j] <- mean(precip[i,j,time=(1:31)])}}，非
→ AndrewShi: 常謝謝andrew大用引導式的方式教我，其實我也比較喜歡
→ AndrewShi: 用這種方式來學習。
→ andrew43: 嗯嗯對
→ andrew43: 等等。第一句不對啊。
→ andrew43: 幹嘛把容器的值在計算前就填入舊值。
→ andrew43: 我的例子不是這麼寫的
推 AndrewShi: 抱歉，andrew大，剛剛才看到你最後的回覆，所以我要改
→ AndrewShi: 成precip1 <- matrix(precip1,1440,400)這樣才對嗎?!
→ andrew43: 創造一個矩陣，裡面先全填成NA或0或某個特定數字都好
→ andrew43: 若填成NA，之後你loop完之後，如果還有NA就表示過程有問
→ andrew43: 題。你先填入舊資料就不能辨視出哪個值可能有問題了。
→ andrew43: 現在你又問是不是先填入precip1，問題是precip1還不存在
→ andrew43: 前不可能，就像a未定義時 a <- a+1 怎麼進行？
推 AndrewShi: 了解XD，所以我的precip1要改成NA_real_，先把這個建
→ AndrewShi: 立出來空的矩陣(precip1)先(隨便)填個值(NA或0)，之後
→ AndrewShi: 再把算出來的平均降雨值填到這個矩陣裡。
→ andrew43: 我是這樣的意思沒錯
推 AndrewShi: andrew大~我想請教你一個問題，如果畫圖的部分我想畫特
→ AndrewShi: 定區域的話(如:東亞)，改image裡的經.緯度範圍和月平
→ AndrewShi: 均降雨裡的經.緯度範圍可以畫的出來，但是之後再疊加世
→ AndrewShi: 界地圖(有設成同樣範圍)的時候畫不出來(它有顯示有多
→ AndrewShi: 少個警告訊息)，但我打warnings( )，它也沒列出警告訊
→ AndrewShi: 息，想請教你這個問題該怎麼解決呢??
→ andrew43: 請不要一直在推文中提出新問題。發文並附code
推 AndrewShi: 好的。