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