【已解决】求助：求最大降水率出现的当地时间

银酱赛高 · 发表于 2013-4-5 21:59:28

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

本帖最后由银酱赛高于 2013-7-18 15:48 编辑

整理了8个时次的grd文件，批量编写了ctl文件，输出没有问题。
资料都是世界时（00，03，06，09，12，15，18，21），要转换成北京时
'reinit'
'open d:\paper\ave_accu.ctl'
'set lon 109 120'
'set lat 20 26'
'set grads off'
'set grid off'
'set cterp on'
'set mpdset cnworld'
'set ylint -1'
'set gxout shaded'
'set parea 1 9.5 0.5 8'
'd:\paper\gs\rgb.gs'
aa=0
while(t <=8)
if(pcp>=aa)
      'c'
      'aa=pcp'
      if(t >=7)
         'd (t -1)*3-16'
      else
         'd (t -1)*3+8'
      endif
endif
t=t +1
endwhile
'd:\paper\gs\cbar_matlab 2 0.7 0'
'printim d:\paper\weixiang.jpg x1000 y800 white'
;
运行后提示  cannot plot color bar:no shading information
但是之前画了几个图都是用一样的shaded参数设置，到底哪一步出问题呢？

mofangbao · 发表于 2013-4-6 00:39:52

楼主把色标之类的额外设置去掉的话，图是否正常，不正常的话是如何不正常，最好图贴出来

银酱赛高 · 发表于 2013-4-6 11:00:52

mofangbao 发表于 2013-4-6 00:39

登录/注册后可看大图

楼主把色标之类的额外设置去掉的话，图是否正常，不正常的话是如何不正常，最好图贴出来

把额外设置都去掉
'reinit'
'open d:\paper\ave_accu.ctl'
'set gxout shaded'
aa=0
while(t <=8)
if(pcp>=aa)
      'c'
      'aa=pcp'
      if(t >=7)
         'd (t -1)*3-16'
      else
         'd (t -1)*3+8'
      endif
endif
t=t +1
endwhile
'printim d:\paper\weixiang.jpg x1000 y800 white'
;
Grads没有报错，但也没有图输出，一片空白

mofangbao · 发表于 2013-4-6 11:14:27

把你程序逻辑和变量含义介绍下吧

银酱赛高 · 发表于 2013-4-6 11:33:04

mofangbao 发表于 2013-4-6 11:14

登录/注册后可看大图

把你程序逻辑和变量含义介绍下吧

ctl批量描述了8个时次（世界时00，03，06，09，12，15，18，21）的平均降水率资料，现在要输出最大降水率出现的当地时间(北京时间)，就是
t                1    2    3    4 5    6    7    8
UTC          00 03 06 09 12 15 18 21
北京时       08 11 14 17 20 23 02 05
当t<7，北京时=(t -1)*3+8；当t>=7，北京时=(t -1)*3-16
所以对t求循环，当pcp大于aa时就清屏，令aa=pcp，输出此时的北京时

mofangbao · 发表于 2013-4-6 12:18:39

把ctl也贴出来吧，这些其实你一开始发帖子就应该贴出来了。
你这个gs有不少问题，aa你定义的是一个脚本变量，只能存放一个单值，而pcp是一个场变量，两种不同的变量不能拿来比较的，如果你pcp是一个只有时间维度的变量，不过没有ctl我看不出来。如果是那样的话，你可以set gxout print,循环set t,然后d变量，结果会保存在result变量中，这时用截取函数把数值取出来，就可以和aa进行比较了，其实你这个用fortran做更好。如果是场变量，我想知道你比较大小的含义是什么，是这个场里面的每一个格点上的值都比某个数大？你的gs实在太矛盾了

银酱赛高 · 发表于 2013-4-6 18:11:23

本帖最后由银酱赛高于 2013-4-6 18:16 编辑

mofangbao 发表于 2013-4-6 12:18

登录/注册后可看大图

把ctl也贴出来吧，这些其实你一开始发帖子就应该贴出来了。
你这个gs有不少问题，aa你定义的是一个脚本变量 ...

没有考虑到变量不同的问题

我的思路是循环8个时次，若降水率pcp>aa就把pcp值赋给aa，并输出此时的时次t，循环最后得出最大降水率及其所在时次，并换算成北京时，但不知道在Grads应该怎么实现~~ grd资料的ctl如下：
dset d:\paper\grd\%h2.grd
options template
undef -9999.9
xdef 45 linear 109          0.25
ydef 25 linear 20          0.25
zdef 1 linear 0          1
tdef 8 linear 00Z01JUN1998  3hr
vars 1
pcp       0 t,y,x precipitation
endvars
我用Fortran尝试了一下，把最大降水率所在时间输出为grd文件，pcp(i,j,t)是降水率，最后求出a(i,j)是最大降水率，t(i,j)是最大降水率所在的时间，最后换算成北京时，Fortran程序如下：

program max_time
implicit none
integer,parameter::M=45
integer,parameter::N=25
integer,parameter::P=8
integer::i,j,k,t(M,N)
integer,parameter::undef=-9999.9
real::pcp(M,N,P),a(M,N)
open(21,file='d:\paper\grd\00.grd',form='binary')
do j=1,25
do i=1,45
read(21) pcp(i,j,1)
end do
end do
close(21)
open(22,file='d:\paper\grd\03.grd',form='binary')
do j=1,25
do i=1,45
read(22) pcp(i,j,2)
end do
end do
close(22)
open(23,file='d:\paper\grd\06.grd',form='binary')
do j=1,25
do i=1,45
read(23) pcp(i,j,3)
end do
end do
close(23)
open(24,file='d:\paper\grd\09.grd',form='binary')
do j=1,25
do i=1,45
read(24) pcp(i,j,4)
end do
end do
close(24)
open(25,file='d:\paper\grd\12.grd',form='binary')
do j=1,25
do i=1,45
read(25) pcp(i,j,5)
end do
end do
close(25)
open(26,file='d:\paper\grd\15.grd',form='binary')
do j=1,25
do i=1,45
read(26) pcp(i,j,6)
end do
end do
close(26)
open(27,file='d:\paper\grd\18.grd',form='binary')
do j=1,25
do i=1,45
read(27) pcp(i,j,7)
end do
end do
close(27)
open(28,file='d:\paper\grd\21.grd',form='binary')
do j=1,25
do i=1,45
read(28) pcp(i,j,8)
end do
end do
close(28)
do j=1,25
do i=1,45
a(i,j)=0.
end do
end do
do j=1,25
do i=1,45
t(i,j)=0
end do
end do
do k=1,8
do j=1,25
do i=1,45
if(pcp(i,j,k)>a(i,j)) then
a(i,j)=pcp(i,j,k)
t(i,j)=k
end if
end do
end do
end do
open(30,file='d:\paper\grd\weixiang.grd',form='binary')
do j=1,25
do i=1,45
if(t(i,j)<7) then
t(i,j)=(t(i,j)-1)*3+8
else
t(i,j)=(t(i,j)-1)*3-16
end if
write(30) t(i,j)
end do
end do
end

复制代码

成功输出后，再对输出文件编写ctl如下：
dset d:\paper\grd\weixiang.grd
undef -9999.9
xdef 45 linear 109          0.25
ydef 25 linear 20          0.25
zdef 1 linear 0          1
tdef 1 linear 00Z01JUN1998  1hr
vars 1
a  0 99  shici
endvars

但Grads一直提示constant field.value=1.4013e-45
然后我就把t(i,j)改成dat输出，都是北京时数据，而且大部分是白天，应该挺正常的；
再输出程序里面的a(i,j),pcp(i,j,k)数组，用Grads画图也是正常的；
在Fortran里面找debug也找不着，实在不知道问题出在哪儿了~

mofangbao · 发表于 2013-4-6 20:21:16

有没有把你的t输出来看看，还有建议a的初始值为-1，把你每一个求t的地方的各个参数输出来看看，思路应该没问题，你哪里的细节弄错了，看看你的公式有没有问题

银酱赛高 · 发表于 2013-4-6 21:43:29

mofangbao 发表于 2013-4-6 20:21

登录/注册后可看大图

有没有把你的t输出来看看，还有建议a的初始值为-1，把你每一个求t的地方的各个参数输出来看看，思路应该没问 ...

有在95行和107行输出过t(i,j)的txt文件来看（一个是时次，一个是北京时），两者关系是对应的；但在相同行位置输出的grd都不能用Grads来打开。于是我编了个Fortran程序把这两个grd再输出为txt来看，结果都是一堆 7.0064923E-45这样的数字。是Fortran程序输出部分的问题吗？

宁小才 · 发表于 2013-4-29 11:52:07

很给力

		自动登录	找回密码
密码			立即注册

[分享资料] 【已解决】求助：求最大降水率出现的当地时间

登录后查看更多精彩内容~