# 已先自 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/cn.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