关于站点数据多时次处理和插值保存问题

Sylvia

       本人之前在论坛获得许多帮助，气象方面的东西说复杂不复杂，说简单也不简单，很多时候就是一些小问题就把大家困住了，于是论坛成了我们求助的对象。我在这里学到很多，有时可以基本解决问题，有时可以得到一些思路和想法。于是，我也想把自己获得的一些方法发一下贴，让一些和我一样的同学们、朋友们能够得到启发和帮助，大家共同学习。小小经验，不足挂齿，希望各位牛人和各路大神多多提点。
      论坛里已经多次提到处理站点数据的方法，但还是有很多人提问关于多时次站点的提取问题，以及有人想要将站点数据保存下来以便之后处理，当然，你们遇到的问题，我和一些人也遇到过，尝试过其实很简单，在此简单做个共享。在此，我要多谢清风的帖子，他关于站点插值过程写的真的很全很好，大家还有看不懂的可以去看这篇帖子，我在此只作一点多时次和插值后的数据保存的补充。多谢。http://bbs.06climate.com/forum.php?mod=viewthread&tid=7357还有陌版主的帖子也很好
http://bbs.06climate.com/forum.php?mod=viewthread&tid=7340，大家都可以借鉴下。
    下面附上处理站点资料多时次处理的fortran程序代码。
      program stngrid
real a3(160,62),a4(160,62),a5(160,62),lat(160),lon(160),x(160,62)
character*8 id(160)
real a
open(3,file='t1606.txt')
  open(4,file='t1607.txt')
open(5,file='t1608.txt')ccccccccccccccccccccc6,7,8月的资料
open(12,file='160.txt')cccccccccccccccccccccccc台站数据
open(11,file='ave.txt')
open(13,file='ave.grd',form='binary')
   
ccccccccccccccc 读数据（经纬度、160站温度）
      read(3,*)((a3(i,j),i=1,160),j=1,62)
      read(4,*)((a4(i,j),i=1,160),j=1,62)
      read(5,*)((a5(i,j),i=1,160),j=1,62)
cccccccccccccccccccccccccc求平均，不求平均，省去即可
ave=0.
ano=0.
      do j=1,62
      do i=1,160
           ave(i)=ave(i)+((a3(i,j)+a4(i,j)+ a5(i,j))/3)/62
      enddo;enddo

  cccccccccc读台站数据   
do i=1,160
read(12,*)a,lat(i),lon(i)
enddo
      write(*,*)' ok 1'
cccccccccccccc写站点数据
      do i=1,62cccccccccccccccccccccccccccc多时次读取，这里是62年，自己改值
do j=1,160cccccccccccccccccccccccccccc160站，使用其他站的资料，也可以照样改值即可
id(j)=char(j)
tim=0.0
nlev=1
nflag=1
write(9)id(j),lat(j),lon(j),tim,nlev,nflag,x(j,i)ccccccccccccccccccccccccccc单时次为x(i)即可
enddoccccccccccccccccccccccccccccccccccccccccc
tim=0.0
nlev=0
nflag=1
      write(9)id(j-1),lat(j-1),lon(j-1),tim,nlev,nflag
      end doccccccccccccccccccccccccccccccccccccccccccc若不是多时次，则此do循环,enddo则不用了
end


因为是多时次，所以在编写格点文件grid.grd和grid.ctl时，对应的时次一定要和站点数据的t时次一致，否则站点映射文件会出错。相对应的grid.grd的fortran程序代码如下：
program  grid
parameter(nx=71,ny=41,nt=62)
real lat(ny),lon(nx)
real s(nx,ny,nt)
open(1,file='grid.grd',form='binary')
lat(1)=15.0
lon(1)=70.0
do j=1,ny-1
   lat(j+1)=lat(j)+1.0
enddo
do i=1,nx-1
   lon(i+1)=lon(i)+1.0
enddo
do k=1,62
  do i=1,nx
    do j=1,ny
      s(i,j,k)=1
enddo
  enddo
enddo
write(1) s
end


相应的ctl分别是station.ctl
dset station.grd
dtype station
stnmap station.map
undef -999.0
title temeprature
tdef 62 1 linear jan1951 1 mo
vars 1
t 0 99 data
endvars


grid.ctl
dset grid.grd
undef -9.99e33
title Sample Grid Data
xdef 71 linear 70 1
ydef 41 linear 15 1
zdef 1 linear 1000 1
tdef 62 linear JAN1951 1yr
vars 1
g 0 99 grid data
endvars


插值画图
'open grid.ctl'
'open station.ctl'
'enable print ave.gmf'
'set vpage 2 10 1 8'
'set cthick 5'
'set lat 15 60'
'set lon 70 150'
'set mpdset hirs cnworld'
i=1
while(i<=62)
'set t 'i''
'define a=oacres(g,t.2,1.5)'
'define a1=maskout(a,g-0.5)'
'define aa=smth9(a1)'
'set lat 16 56'
'set lon 75 135'
'set xlopts 1 4 0.12'
'set ylopts 1 5 0.12'
'set gxout contour'
'set clab forced'
'set ccolor rainbow'
'set grads off'
'set grid off'
'set cint 0.3'
'cnbasemap aa*0.1'
'run china -p'
'print'
'c'
i=i+1
endwhile
'disable print'
；
插值函数和图形美化的一些方法在此就不多说了，论坛帖子里都能找得到。
图形出来了，于是有人不想出图，而是把图保存下来，于是我们用到了set fwrite 这个命令。之前我发过一篇帖子，我使用了这个命令但是解决不了问题。在清风的帖子里也没有。有人会用兰溪之水的帖子gradsascii函数转成txt格式的，但是想直接放到fortran里处理还要转换比较麻烦，当然，清风和兰溪之水的帖子介绍的很详细，推荐大家去参考，我也是借用他们的帖子得出的方法，呵呵。后来我找了其他的资料，尝试了几次之后，发现这个命令还是可以用的，发现set fwrite 这个命令还是很强大的。这里只是简单的调整下set fwrite 的顺序即可。
'open grid.ctl'
'open station.ctl'
'set lat 15 55'
'set lon 70 140'
i=1
while (i<=62)
'set t 'i''
'set gxout fwrite'cccccccccccccccccccccccccccccc一般set fwrite 和set gxout fwrite 是反过来写的，在这里在颠倒顺序，并且在多次里需要将set fwrite写到时次的循环里，才能得到准确的62个时次的数据，否则只有同样的值或者一个时次的值。
'set fwrite ave.grd'

'define a=oacres(g,t.2,1.5)'
'define a1=maskout(a,g-0.5)'
'define aa=smth9(a1)'
'd aa*0.1'
i=i+1
endwhile
'disable fwrite'
;
  就不多说废话了，到此为止。。。。后面可能会将处理其他站点数据的方法再放上来还有一些处理图形的技巧，都是在大家的帮助下边学边收获的经验，希望大家共同分享。

兰溪之水

grid数据可以只写一个时次，调用时用第一时次数据就行~关于fwrite个人觉得下面写就可以，你说不可以，你能再试试么？
'open grid.ctl'
'open station.ctl'
'set lat 15 55'
'set lon 70 140'

'set fwrite ave.grd'
'set gxout fwrite'
i=1
while (i<=62)
'set t 'i''
'define a=oacres(g(t=1),t.2,1.5)'
'define a1=maskout(a,g(t=1)-0.5)'
'define aa=smth9(a1)'
'd aa*0.1'
i=i+1
endwhile
'disable fwrite'
;

mofangbao

改天实验一下，很感谢楼主的分享，果然说话很算数，哈哈

Sylvia

mofangbao 发表于 2012-9-16 17:03

                               
登录/注册后可看大图

改天实验一下，很感谢楼主的分享，果然说话很算数，哈哈

我会再接再厉哈，还有很多问题，一个一个解决，向你们看齐。。。。。

做个霸气的木头

真是辛苦楼主了，这么无私的分享给大家~

Sylvia

兰溪之水 发表于 2012-9-16 22:30

                               
登录/注册后可看大图

grid数据可以只写一个时次，调用时用第一时次数据就行~关于fwrite个人觉得下面写就可以，你说不可以，你能再 ...

我之前试过，是不可以的，但是这也许是每个人电脑和grads的不同造成的。之前也遇到过这样的问题，可能在别的电脑上能行。也可能过几天这么试就可以了。至于那个时次的问题，你那样调用第一个时次是可以的，我也看到过，但是我尝试过，我的不行，可能是数据的问题。李老师说过时次最好一定要一致，以免出错，还是保险点好。

冰点

你好！你看编写格点文件简写成这样可好？
program  griddata
integer,parameter::nx=71,ny=41,nt=62
real grid
integer i,j,k
open(11,file='e:\grid.grd',form='binary')
grid=1.0
do k=1,nt
    do i=1,nx
        do j=1,ny
            write(11)grid
        enddo
    enddo
enddo
close(11)
stop
end

Sylvia

冰点 发表于 2012-9-19 09:46

                               
登录/注册后可看大图

你好！你看编写格点文件简写成这样可好？
program  griddata
integer,parameter::nx=71,ny=41,nt=62

这样子是默认经纬度吗？没有设置精度，也没有设置起止经纬度呀

mofangbao

Sylvia 发表于 2012-9-20 09:37

                               
登录/注册后可看大图

这样子是默认经纬度吗？没有设置精度，也没有设置起止经纬度呀

经纬度啥的是在ctl中设置的啦，对于格点数据，fortran只要写入数据就行啦，具体数据代表的范围由ctl确定

Sylvia

mofangbao 发表于 2012-9-20 09:43

                               
登录/注册后可看大图

经纬度啥的是在ctl中设置的啦，对于格点数据，fortran只要写入数据就行啦，具体数据代表的范围由ctl确定

o ,ctl这么强大呀，呵呵，好像是的呢。谢谢哈