# 已先自 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檔 降雨

非常感謝andrew大，願意花這麼多時間幫我解答，我今天也會好好研究程式碼，如果對於你寫的程式碼有疑惑的話方便能再請教你嗎XD??abdrew大~我想請教apply裡的c(1,2)是指把每列.每行(每個經度.緯度)的降雨值都帶到function裡面的意思對嗎?!另外想請問你是怎麼把時間(月)放到precip.ave.monthly的第一個維度裡的呢??

這應該是因為tapply()一次有12個值，所以apply自動把它安排在第一個維度。沒有哪段碼刻意要求達到這種結果。再接一個aperm(..., c(2,3,1))轉換過來就好了另外，單用apply算降雨年平均後只有二維度，並沒有時間.

了解，真的非常感謝andrew大，讓我學到很多~andrew大~想再請教你如果想用迴圈的概念來寫的話，我想到的寫法是:for(i in (1:1440)){for(j in (1:400)){mean(precip[i,j,time=(1:31)])}}但跑出來也只有一個值，所以想請教你我的迴圈寫法(概念)是哪裡有出錯嗎??

這樣是各別算出每個地點一月份的平均雨量，一次算一個。那還不如寫成apply(precip[,,1:31], c(1,2), mean)你若想用for loop一次算一個地點，要把結果填到一個容器中，例如一個400*1440矩陣中。

了解，非常感謝~

嗯嗯對等等。第一句不對啊。幹嘛把容器的值在計算前就填入舊值。我的例子不是這麼寫的

抱歉，andrew大，剛剛才看到你最後的回覆，所以我要改成precip1 <- matrix(precip1,1440,400)這樣才對嗎?!

創造一個矩陣，裡面先全填成NA或0或某個特定數字都好若填成NA，之後你loop完之後，如果還有NA就表示過程有問題。你先填入舊資料就不能辨視出哪個值可能有問題了。現在你又問是不是先填入precip1，問題是precip1還不存在前不可能，就像a未定義時 a <- a+1 怎麼進行？

了解XD，所以我的precip1要改成NA_real_，先把這個建立出來空的矩陣(precip1)先(隨便)填個值(NA或0)，之後再把算出來的平均降雨值填到這個矩陣裡。

我是這樣的意思沒錯

andrew大~我想請教你一個問題，如果畫圖的部分我想畫特定區域的話(如:東亞)，改image裡的經.緯度範圍和月平均降雨裡的經.緯度範圍可以畫的出來，但是之後再疊加世界地圖(有設成同樣範圍)的時候畫不出來(它有顯示有多少個警告訊息)，但我打warnings( )，它也沒列出警告訊息，想請教你這個問題該怎麼解決呢??

請不要一直在推文中提出新問題。發文並附code

好的。