# 已先自 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://webptt.com/m.aspx?n=bbs/R_Language/M.1540258357.A.42A.html precip 是一個三維陣列，維度按順序分別是經度、緯度和時間。 apply(precip, c(1, 2), mean) 的意思是， 對 precip 做平均，方式是同經度同緯度而不同時間的值取成一組資料的平均。 這裡的 c(1, 2) 就是指第一和第二個維度，即同經度同緯度的意思。 因此，其結果就由365「層」平均成 1 層。 上面如果可以理解，那麼 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 不對。 你似乎不太熟悉一些矩陣或陣列的操作。 請研究以下例子。 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