[问题] 循环问题(ncdf档)

楼主: AndrewShi (没有妳的我)   2019-05-09 18:52:14
[问题类型]:
程式咨询(我想用R 做某件事情,但是我不知道要怎么用R 写出来)
[软件熟悉度]:
入门(写过其他程式,只是对语法不熟悉)
[问题叙述]:
各位大大好,小弟我目前有2个矩阵(原始为list)的数值资料,其中一个为时间的
index(第几笔),另一个为各别分配的权重,我要做的事是将这2个矩阵的资讯套用在一笔
降雨资料(netcdf档)上,由下图的2个矩阵(时间index.权重)的资讯为例,也就是我要将
原始的降雨资料的第一天(笔)以第1天的降雨值(每个网格点)x0.8334249,加上第9830天的
降雨值x0.12252973来取代,以此类推。
而小弟我目前卡在不知道该怎么把矩阵对应的关系(第1天的第1个时间index乘上第一
个权重)放入循环中,下面的程式码中是以analog.indices表示时间的index,weights表
示权重,可能程式码的逻辑非常怪异,我上网查或许用apply系列的指令较适合,但我也
较少用apply系列的指令,因此较不熟悉,还烦请大大们指点和较详细的说明,也非常欢
迎引导式教学,谢谢。
2个矩阵以及降雨的资料放在此:http://0rz.tw/JI056
https://imgur.com/Q4rRKi9 (时间index)
https://imgur.com/YJFXtr5 (权重)
[程式范例]:
library(ncdf4)
library(data.table)
memory.limit(size=50000)
analog.indices <- A[[1]]
analog.indices <- matrix(unlist(analog.indices), nrow=length(analog.indices),
byrow=T)
weights <- A[[2]]
weights <- matrix(unlist(weights), nrow=length(weights), byrow=T)
HIRAM_WRF_data <- nc_open("C:\\Users\\TOM\\Desktop\\R(数据库)\\WRF(动力降尺度
资料)\\T2WHIRAM_c384_amip\\197901-200512_pr_axis_time_domain.nc")
print(HIRAM_WRF_data)
hiram_wrf_lon <- ncvar_get(HIRAM_WRF_data,"lon")
hiram_wrf_lat <- ncvar_get(HIRAM_WRF_data,"lat")
hiram_wrf_time <- ncvar_get(HIRAM_WRF_data,"time")
hiram_wrf_pr <- ncvar_get(HIRAM_WRF_data,"pr")
pr <- array(NA,dim=c(length(hiram_wrf_lon),length(hiram_wrf_lat),2))
for(i in analog.indices[1,i]){
for(w in weights[1,w]){
if(i==w){
pr[,,1:2] <-
ncvar_get(HIRAM_WRF_data,"pr",start=c(1,1,i),count=c(41,77,1))*w
}
}
}
View(pr[,,1])
[环境叙述]:
[关键字]:
循环 ncdf档
作者: andrew43 (讨厌有好心推文后删文者)   2019-05-09 19:08:00
analog.indices的V1有重复,对吗?例如V1==10另外,这用for会好写很多。晚点再看看权重表格如何和其他对应?说清楚一些
楼主: AndrewShi (没有妳的我)   2019-05-09 20:41:00
andrew大~对,但应该不多,analog.indices那个表格最左边那一列是原始资料的天数,右边这两列则是所有资料里和这天最相近的2天。权重的表格是和analog.indices的表格相对应,也就是两者的[1,V1]做相乘,两者的[1,V2]做相乘,再相加,只是做相乘的时候analog.indices的[1,V1]不是取那个数字,而是取那天的降雨资料出来和权重的[1,V1]这个值做相乘,[1,V2]也是一样。
作者: andrew43 (讨厌有好心推文后删文者)   2019-05-09 21:28:00
那么,遇到第10天时怎么算?能把实际数字写出来吗?
楼主: AndrewShi (没有妳的我)   2019-05-09 22:15:00
一样呀,第10天就是拿第10天的降雨值乘上权重[10,V1](0.764485),加上第4407天的降雨值乘上权重[10,V2](0.203293)相加而得。
作者: andrew43 (讨厌有好心推文后删文者)   2019-05-09 23:31:00
但有2个第10天不是吗?有重复要怎么处理?
楼主: AndrewShi (没有妳的我)   2019-05-10 01:20:00
那是第10天被选了2次,意思是原始资料的第10天和第11天很像,所以原始资料在第10天选了自己本身,在第11天时也选了第10天,右边这两列的数字是选出跟原本资料(最左列)那天最相近的2天出来,所以基本上会挑出自己本身(第一列和第二列数字几乎一样);然而原始资料的第10天(最左边的列)是没有重复的。
作者: andrew43 (讨厌有好心推文后删文者)   2019-05-10 08:38:00
ok。那么,雨量资料有9862天,但对应表格有10227列,而你上述对应表格的列号就是对应第几天。为什么?例如,对应表格第10227列的作用是什么?或是你新创出的雨量资料要有10227层?是的话,参考 https://ideone.com/2S9KgX
楼主: AndrewShi (没有妳的我)   2019-05-10 11:11:00
表格的列号正常的话也只有9862列,后面是我新创的,忘记删掉,抱歉造成你理解上的误会,再次感谢andrew大,code的部分我会再好好的研究。
作者: andrew43 (讨厌有好心推文后删文者)   2019-05-10 11:18:00
那把新array的dim改一下就可以了
楼主: AndrewShi (没有妳的我)   2019-05-10 13:20:00
恩恩,andrew大,想再请教你如果我每天的analog.indices和weights都有30个的话,col.names和for循环的部分除了自己补足剩下的,该怎么修改循环中$i1.$i2.$w1.$w2的部分让它自动加总30天呢??
作者: andrew43 (讨厌有好心推文后删文者)   2019-05-10 13:36:00
最直接的想法是双层循环,外面那层照旧走直的,里面那层走横的。另一种可以是把权重表格重新整理成特定顺序,再利用矩阵乘法做加权。速度应该会快很多但较不直观。
楼主: AndrewShi (没有妳的我)   2019-05-10 15:16:00
我想到的是第一种(双层循环),不知道改成这样对不对,还请andrew大指点。https://imgur.com/ci54KFs
作者: andrew43 (讨厌有好心推文后删文者)   2019-05-10 15:19:00
不对喔。错在w那层没有用处,变成只是同一件事做很多次.晚点我写给你看看
楼主: AndrewShi (没有妳的我)   2019-05-10 15:19:00
因为不知道怎么把$i1.$i2.$w1.$w2写成循环,所以把i1.i2.w1.w2全改成i和w。
作者: andrew43 (讨厌有好心推文后删文者)   2019-05-10 16:19:00
楼主: AndrewShi (没有妳的我)   2019-05-10 18:04:00
原来是要改成这样,非常感谢andrew大,最后有2个小问题想请教你,一个是.[1:9862]最前面的那个.是代表前面读进来的那个csv档吗?!另一个是col.names如果i有1~30的话是只能一个一个打(命名)吗??
作者: andrew43 (讨厌有好心推文后删文者)   2019-05-10 18:06:00
1. 是2. 此例命名不必要,随便取成as.character(1:30)也行
楼主: AndrewShi (没有妳的我)   2019-05-10 23:38:00
了解,再次感谢andrew大的解答和每次的帮助~

Links booklink

Contact Us: admin [ a t ] ucptt.com