PTT數位生活區 / R_Language

Re: [問題] 月平均資料

看板R_Language作者 (討厭有好心推文後刪文者)時間6年前 (2018/10/23 09:32), 6年前編輯推噓13(13推 0噓 54→)
留言67則, 2人參與, 6年前最新討論串2/2 (看更多)
# 已先自 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
推
10/23 11:32, 6年前 , 1F
非常感謝andrew大，願意花這麼多時間幫我解答，我今天
10/23 11:32, 1F
→
10/23 11:32, 6年前 , 2F
也會好好研究程式碼，如果對於你寫的程式碼有疑惑的話
10/23 11:32, 2F
→
10/23 11:32, 6年前 , 3F
方便能再請教你嗎XD??
10/23 11:32, 3F
推
10/23 15:24, 6年前 , 4F
abdrew大~我想請教apply裡的c(1,2)是指把每列.每行(每
10/23 15:24, 4F
→
10/23 15:24, 6年前 , 5F
個經度.緯度)的降雨值都帶到function裡面的意思對嗎?!
10/23 15:24, 5F
precip 是一個三維陣列，維度按順序分別是經度、緯度和時間。 apply(precip, c(1, 2), mean) 的意思是， 對 precip 做平均，方式是同經度同緯度而不同時間的值取成一組資料的平均。 這裡的 c(1, 2) 就是指第一和第二個維度，即同經度同緯度的意思。 因此，其結果就由365「層」平均成 1 層。
→
10/23 15:24, 6年前 , 6F
另外想請問你是怎麼把時間(月)放到precip.ave.monthly
10/23 15:24, 6F
→
10/23 15:24, 6年前 , 7F
的第一個維度裡的呢??
10/23 15:24, 7F
上面如果可以理解，那麼 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
推
10/24 14:47, 6年前 , 8F
感謝andrew大詳細的解答，上面的敘述我能理解，不過我
10/24 14:47, 8F
→
10/24 14:47, 6年前 , 9F
很好奇為什麼在單用apply算降雨年平均的時候時間是在
10/24 14:47, 9F
→
10/24 14:47, 6年前 , 10F
降雨的第三個維度，而在算月平均的時候，時間變成降雨(
10/24 14:47, 10F
→
10/24 14:47, 6年前 , 11F
precip.ave.monthly)的第一個維度，而經.緯度則變成第
10/24 14:47, 11F
→
10/24 14:47, 6年前 , 12F
二.三個維度，我看不出來哪一段程式碼是在做這樣的處理
10/24 14:47, 12F
→
10/24 14:47, 6年前 , 13F
。
10/24 14:47, 13F
→
10/24 16:59, 6年前 , 14F
這應該是因為tapply()一次有12個值，所以apply自動把它
10/24 16:59, 14F
→
10/24 17:01, 6年前 , 15F
安排在第一個維度。沒有哪段碼刻意要求達到這種結果。
10/24 17:01, 15F
→
10/24 17:03, 6年前 , 16F
再接一個aperm(..., c(2,3,1))轉換過來就好了
10/24 17:03, 16F
→
10/24 17:05, 6年前 , 17F
另外，單用apply算降雨年平均後只有二維度，並沒有時間.
10/24 17:05, 17F
推
10/24 22:25, 6年前 , 18F
了解，真的非常感謝andrew大，讓我學到很多~
10/24 22:25, 18F
推
10/25 16:22, 6年前 , 19F
andrew大~想再請教你如果想用迴圈的概念來寫的話，我想
10/25 16:22, 19F
→
10/25 16:22, 6年前 , 20F
到的寫法是:for(i in (1:1440)){
10/25 16:22, 20F
→
10/25 16:22, 6年前 , 21F
for(j in (1:400)){
10/25 16:22, 21F
→
10/25 16:22, 6年前 , 22F
mean(precip[i,j,time=(1:31)])}}
10/25 16:22, 22F
→
10/25 16:22, 6年前 , 23F
但跑出來也只有一個值，所以想請教你我的迴圈寫法(概
10/25 16:22, 23F
→
10/25 16:22, 6年前 , 24F
念)是哪裡有出錯嗎??
10/25 16:22, 24F
→
10/25 16:47, 6年前 , 25F
這樣是各別算出每個地點一月份的平均雨量，一次算一個。
10/25 16:47, 25F
→
10/25 16:48, 6年前 , 26F
那還不如寫成apply(precip[,,1:31], c(1,2), mean)
10/25 16:48, 26F
→
10/25 16:52, 6年前 , 27F
你若想用for loop一次算一個地點，要把結果填到一個容器
10/25 16:52, 27F
→
10/25 16:53, 6年前 , 28F
中，例如一個400*1440矩陣中。
10/25 16:53, 28F
推
10/25 17:09, 6年前 , 29F
了解，非常感謝~
10/25 17:09, 29F
推
10/26 17:28, 6年前 , 30F
andrew大~不好意思，想再請教你一個基本問題，把結果
10/26 17:28, 30F
→
10/26 17:28, 6年前 , 31F
填到一個矩陣裡，matrix(mean(precip[i,j,time=(1:31),
10/26 17:28, 31F
→
10/26 17:28, 6年前 , 32F
1440,400))這樣寫對嗎??
10/26 17:28, 32F
不對。 你似乎不太熟悉一些矩陣或陣列的操作。 請研究以下例子。 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
推
10/27 17:21, 6年前 , 33F
我懂了，首先先創一個數字由1~24，2列3行共4個(矩陣)
10/27 17:21, 33F
→
10/27 17:21, 6年前 , 34F
的陣列a，之後再創一個矩陣m，列和行的數目和陣列a的
10/27 17:21, 34F
→
10/27 17:21, 6年前 , 35F
第一和第二個維度一樣(2列3行)，最後再把陣列a每個列和
10/27 17:21, 35F
→
10/27 17:21, 6年前 , 36F
行各別的值(共4個)相加取平均後放到矩陣m裡(如:(1+7+13
10/27 17:21, 36F
→
10/27 17:21, 6年前 , 37F
+19)/4=10)，我這樣的理解應該沒有錯吧 :)?!
10/27 17:21, 37F
推
10/27 17:30, 6年前 , 38F
而我的程式應該改為:precip1 <- matrix(precip,1440,40
10/27 17:30, 38F
→
10/27 17:30, 6年前 , 39F
0)
10/27 17:30, 39F
→
10/27 17:30, 6年前 , 40F
for(i in (1:1440)){
10/27 17:30, 40F
→
10/27 17:30, 6年前 , 41F
for(j in 1:400)){
10/27 17:30, 41F
→
10/27 17:30, 6年前 , 42F
precip1[i,j] <- mean(precip[i,j,time=(1:31)])}}，非
10/27 17:30, 42F
→
10/27 17:30, 6年前 , 43F
常謝謝andrew大用引導式的方式教我，其實我也比較喜歡
10/27 17:30, 43F
→
10/27 17:30, 6年前 , 44F
用這種方式來學習。
10/27 17:30, 44F
→
10/27 18:28, 6年前 , 45F
嗯嗯對
10/27 18:28, 45F
→
10/27 22:45, 6年前 , 46F
等等。第一句不對啊。
10/27 22:45, 46F
→
10/27 22:45, 6年前 , 47F
幹嘛把容器的值在計算前就填入舊值。
10/27 22:45, 47F
→
10/27 22:47, 6年前 , 48F
我的例子不是這麼寫的
10/27 22:47, 48F
推
11/02 16:31, 6年前 , 49F
抱歉，andrew大，剛剛才看到你最後的回覆，所以我要改
11/02 16:31, 49F
→
11/02 16:31, 6年前 , 50F
成precip1 <- matrix(precip1,1440,400)這樣才對嗎?!
11/02 16:31, 50F
→
11/02 16:56, 6年前 , 51F
創造一個矩陣，裡面先全填成NA或0或某個特定數字都好
11/02 16:56, 51F
→
11/02 16:58, 6年前 , 52F
若填成NA，之後你loop完之後，如果還有NA就表示過程有問
11/02 16:58, 52F
→
11/02 16:59, 6年前 , 53F
題。你先填入舊資料就不能辨視出哪個值可能有問題了。
11/02 16:59, 53F
→
11/02 16:59, 6年前 , 54F
現在你又問是不是先填入precip1，問題是precip1還不存在
11/02 16:59, 54F
→
11/02 17:01, 6年前 , 55F
前不可能，就像a未定義時 a <- a+1 怎麼進行？
11/02 17:01, 55F
推
11/02 18:16, 6年前 , 56F
了解XD，所以我的precip1要改成NA_real_，先把這個建
11/02 18:16, 56F
→
11/02 18:16, 6年前 , 57F
立出來空的矩陣(precip1)先(隨便)填個值(NA或0)，之後
11/02 18:16, 57F
→
11/02 18:16, 6年前 , 58F
再把算出來的平均降雨值填到這個矩陣裡。
11/02 18:16, 58F
→
11/02 19:04, 6年前 , 59F
我是這樣的意思沒錯
11/02 19:04, 59F
推
11/07 17:19, 6年前 , 60F
andrew大~我想請教你一個問題，如果畫圖的部分我想畫特
11/07 17:19, 60F
→
11/07 17:19, 6年前 , 61F
定區域的話(如:東亞)，改image裡的經.緯度範圍和月平
11/07 17:19, 61F
→
11/07 17:19, 6年前 , 62F
均降雨裡的經.緯度範圍可以畫的出來，但是之後再疊加世
11/07 17:19, 62F
→
11/07 17:20, 6年前 , 63F
界地圖(有設成同樣範圍)的時候畫不出來(它有顯示有多
11/07 17:20, 63F
→
11/07 17:20, 6年前 , 64F
少個警告訊息)，但我打warnings( )，它也沒列出警告訊
11/07 17:20, 64F
→
11/07 17:20, 6年前 , 65F
息，想請教你這個問題該怎麼解決呢??
11/07 17:20, 65F
→
11/08 10:36, 6年前 , 66F
請不要一直在推文中提出新問題。發文並附code
11/08 10:36, 66F
推
11/08 12:01, 6年前 , 67F
好的。
11/08 12:01, 67F
‣ 更多 andrew43 的文章
文章代碼(AID): #1RpderGg (R_Language)
討論串 (同標題文章)
本文引述了以下文章的的內容：
完整討論串 (本文為第 2 之 2 篇)：
排序： 最舊先 | 最新先 | 留言數
在新視窗開啟完整討論串 (共2篇)
R_Language 近期熱門文章
更多 近期熱門文章 >>
PTT數位生活區 即時熱門文章
更多 即時熱門文章 >>
‣ 更多 andrew43 的文章
文章代碼(AID): #1RpderGg (R_Language)