DRADS插值时遇到warning from OACRES:less than two station.求指教！

风之叶 · 发表于 2014-4-30 14:48:14

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

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

x

我在用GRADS做站点数据插值到格点时，运行总是出现如题的错误提示。求大家帮忙指导一下！
站点数据描述的ctl文件：数据时间从1958年1月到2009年12月
dset d:\hgt100_00.grd
dtype station
stnmap d:\hgt.map
undef -999.0
title hgt
tdef 624 linear Jan1958 1mo
vars 1
h 0 99 hgt data
endvars

格点文件生成的FORTRAN程序：
parameter(nx=71,ny=41,nt=624)
real lat(ny),lon(nx)
real g(nx,ny,nt)
open(1,file='d:\grid.grd',form='binary')
lat(1)=15.0
lon(1)=70.0
do j=1,ny-1
lat(j+1)=lat(j)+1.0
end do
do i=1,nx-1
lon(i+1)=lon(i)+1.0
end do
do t=1,624
do i=1,nx
  do j=1,ny
  g(i,j,t)=1
  end do
  end do
end do
do t=1,nt
  write(1) ((g(i,j,t),i=1,nx),j=1,ny)
  end do
  end
格点文件相对应的数据描述文件：
dset d:\grid.grd
undef -999.0
title Sample GRID Data
xdef 71 linear 70 1
ydef 41 linear 15 1
zdef 1 linear 100 1
tdef 624 linear Jan1958 1mo
vars 1
g 0 99 grid data
endvars

插值的gs文件：
'open d:\grid.ctl'
'open d:\hgt.ctl'
'set grads off'
'enable print d:\hgt100.gmf'
'set lon 70 140'
'set lat 15 55'
'set mpdset cn cnriver'
'define a=oacres(g,h.2)'
'define a1=maskout(a,g-0.5)'
'define aa=smth9(a1)'
'set xlopts 1 10 0.18'
'set ylopts 1 10 0.18'
'set gxout shaded'
'run d:\rgbset.gs'
'set clevs 0 10 50 100 200 300 400 500 600'
'set ccolor rainbow'
'd aa'
'cbarn.gs'
'set cthick 8'
'set clopts 16 0.1'
'set gxout contour'
'set cint 50'
'd aa'
pull dummy
river(15,4,4)
'run d:\southsea.gs'
'printim d:\hgt100.png white'
' print'
'disable print'
;
求各位指导一下，谢谢！

lqouc · 发表于 2014-4-30 15:54:06

贴一下你站点数据的f文件啦~

风之叶 · 发表于 2014-4-30 18:20:52

lqouc 发表于 2014-4-30 15:54
贴一下你站点数据的f文件啦~

f:\QQ截图20140430181625 这是我拿到手的站点数据，后来处理成二进制了。顺便问你一下，什么是f文件哈？

风之叶 · 发表于 2014-4-30 18:22:18

http://bbs.06climate.com/forum.php?mod=attachment&aid=Mjc1NDB8Y2Y1ZjY2MjVkZGU5YmFmZmY3MWJlMzg5MGM1OGVjNmZ8MTc1MDQ5NDg0Mw%3D%3D&request=yes&_f=.png

风之叶 · 发表于 2014-4-30 19:45:39

lqouc 发表于 2014-4-30 15:54
贴一下你站点数据的f文件啦~

main.f90
program main

use stations ! 使用module子模块，该模块必须与这个主程序放在同一个project下。

implicit none

integer unfyrstart,unfyrend,rlyrstart,rlyrend,nlev,nflag
character*1500::stationline=''
character*9 tempstid
real tim
character*8 stid(150)
character*3 stationnums
integer point,isearch
integer stationnum
integer i,j,k
integer yr,mo,iyr,imo
integer tempa,tempb,tempc,tempd,ioa,iob
real pmeta,pmetb,pmetc,pmetd
real,allocatable::dat(:,:,:),validat(:,:,:)
call filling_coordinate

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 需要修改的地方 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
open(11,file='d:\hgt\00z\hgt100_00.txt',status='old')
!  更改为需要读取的文件路径和名称
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 需要修改的地方 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

open(12,file='Station_Info_in_this_file.txt')    ! 从数据文件中读取台站号后生成的台站信息文件，默认和程序文件在同一路径。

! 判断数据中包含的总台站数
read(11,"(1500a)") stationline
stationnum=1
point=8
tempstid=stationline(point+1:point+9)
do while(tempstid(9:9)/='')
stid(stationnum)=tempstid(5:9)
stationnum=stationnum+1
point=point+9
tempstid=stationline(point+1:point+9)
end do
stationnum=stationnum-1
call num2char(stationnums,stationnum)
write(*,"(a,i3)") "这个数据文件中的台站数为：",stationnum

! 输出数据文件中出现的所有台站信息列表
write(12,"(a)") "序号站号          站名                   经度             纬度          "
do i=1,stationnum
do isearch=1,150
if(stid(i)==coordinate(isearch)%stid) then
exit
end if
end do
write(12,"(i3,4x, a8,11x,a,f20.5,f20.5)")i, stid(i),coordinate(isearch)%name,coordinate(isearch)%lambda,coordinate(isearch)%fi
end do
write(*,"(3a)") "台站信息生成并写入完毕！"

! 自动判读数据文件中的统一起止时间和实际起止时间，与"Data Details.xls"中的信息一致。
read(11,"(1x,i4,5x,f7.1)")tempa,pmeta
unfyrstart=tempa
if(pmeta /=99999.0)then
rlyrstart=unfyrstart
else
do
read(11,"(1x,i4,5x,f7.1)")tempb,pmetb
if(pmetb /=99999.0)then
rlyrstart=tempb
goto 101
endif
enddo
endif

101 do
read(11,"(1x,i4,5x,f7.1)",iostat=ioa)tempc,pmetc
if(ioa /=0)then
unfyrend=tempc
rlyrend=tempc
goto 103
else
if(pmetc ==99999.0)then
rlyrend=tempc-1
goto 102
endif
endif
enddo
102 do
read(11,"(1x,i4,5x,f7.1)",iostat=iob)tempd,pmetd
if(iob /=0)then
unfyrend=tempd
goto 103
endif
enddo
103 continue
write(*,"(a, i4,a,i4)")"本数据文件的统一起止时间为：", unfyrstart, " - ", unfyrend
write(*,"(a, i4,a,i4)")"          实际起止时间为：", rlyrstart, " - ", rlyrend

!  将数据读入数组
rewind(11)
read(11,*)
allocate(dat(stationnum,unfyrstart:unfyrend,12))
allocate(validat(stationnum,rlyrstart:rlyrend,12))
do iyr=unfyrstart,unfyrend
do imo=1,12
read(11,"(1x,i4,1x,i2.2,2x,"//stationnums//"f9.1)") yr,mo,dat(:,iyr,imo)
end do
end do
validat(1:stationnum,rlyrstart:rlyrend,1:12)=dat(1:stationnum,rlyrstart:rlyrend,1:12)
write(*,"(a)") "所有数据及有效数据已分别读入数据 'dat' 和 'validat' "
write(*,"(a,i3,a,i4,a,i4,a)") "  数组'dat' 的维数（台站数，统一起止年份，月份）为：(1:", stationnum,", ",unfyrstart,":",unfyrend,", 1:12)"
write(*,"(a,i3,a,i4,a,i4,a)") "  数组'validat' 的维数（台站数，实际起止年份，月份）为：(1:", stationnum,", ",rlyrstart,":",rlyrend,", 1:12)"

   OPEN(15,FILE='d:\hgt100_00.grd',form='binary')

do j=rlyrstart,rlyrend
do k=1,12
tim=0.0
nlev=1
nflag=1
do i=1,stationnum
do isearch=1,150
if(stid(i)==coordinate(isearch)%stid) then
exit
end if
end do
write(15) STid(i),coordinate(isearch)%lambda,coordinate(isearch)%fi,tim,nlev,nflag,validat(i,j,k)
enddo
nlev=0
write(15) STid(i-1),coordinate(isearch-1)%lambda,coordinate(isearch-1)%fi,tim,nlev,nflag
enddo
enddo

close(15)

!!!!!!!!!接下来，就可以用读入的数据进行分析运算了!!!!!!!!!!!!!!!!!!!!!!!!

end program

subroutine num2char(goal,num)
implicit none
character*3 goal
character*1 ab(3)
integer num,a,b,c
a=num/100
b=num/10-10*a
c=num-10*b-100*a
ab(1)=char(a+48)
ab(2)=char(b+48)
ab(3)=char(c+48)
goal=ab(1)//ab(2)//ab(3)
end subroutine

station.f90
module stations

implicit none
type station
      character*5::stid
      real::lambda
      real::fi
      character*10 name
end type
type(station)::coordinate(150)=station("    ",0.0,0.0,' ')
contains

subroutine filling_coordinate
      implicit none
      character*1 bullshit(100)
      real::lambda_hour=0.0
      real::lambda_minute=0.0
      real::lambda_second=0.0
      real::fi_hour=0.0
      real::fi_minute=0.0
      real::fi_second=0.0
      integer stindex

      open(11,file='station_coordinate.dat',status='old')
      do stindex=1,121
read(11,"(80a)") bullshit(1:80)
coordinate(stindex)%stid=bullshit(1)//bullshit(2)//bullshit(3)//bullshit(4)//bullshit(5)
call char2num__3(bullshit(7:9),lambda_hour)
call char2num__2(bullshit(12:13),lambda_minute)
call char2num__2(bullshit(15:16),lambda_second)
coordinate(stindex)%lambda=sixty2ten(lambda_hour,lambda_minute,lambda_second)
call char2num__2(bullshit(19:20),fi_hour)
call char2num__2(bullshit(23:24),fi_minute)
call char2num__2(bullshit(26:27),fi_second)
coordinate(stindex)%fi=sixty2ten(fi_hour,fi_minute,fi_second)
      end do

      coordinate(1)%name='海拉尔'
      coordinate(2)%name='嫩江'
      coordinate(3)%name='齐齐哈尔'
      coordinate(4)%name='伊春'
      coordinate(5)%name='索伦'
      coordinate(6)%name='哈尔滨'
      coordinate(7)%name='阿勒泰'
      coordinate(8)%name='塔城'
      coordinate(9)%name='克拉玛依'
      coordinate(10)%name='北塔山'
      coordinate(11)%name='伊宁'
      coordinate(12)%name='乌鲁木齐'
      coordinate(13)%name='阿克苏'
      coordinate(14)%name='库车'
      coordinate(15)%name='库尔勒'
      coordinate(16)%name='喀什'
      coordinate(17)%name='若羌'
      coordinate(18)%name='和田'
      coordinate(19)%name='民丰'
      coordinate(20)%name='茫崖'
      coordinate(21)%name='哈密'
      coordinate(22)%name='额济纳旗'
      coordinate(23)%name='马崇山'
      coordinate(24)%name='敦煌'
      coordinate(25)%name='巴彦诺尔公'
      coordinate(26)%name='酒泉'
      coordinate(27)%name='张掖'
      coordinate(28)%name='民勤'
      coordinate(29)%name='格尔木'
      coordinate(30)%name='都兰'
      coordinate(31)%name='西宁'
      coordinate(32)%name='榆中'
      coordinate(33)%name='二连浩特'
      coordinate(34)%name='海流图'
      coordinate(35)%name='呼和浩特'
      coordinate(36)%name='临河'
      coordinate(37)%name='东胜'
      coordinate(38)%name='银川'
      coordinate(39)%name='太原'
      coordinate(40)%name='邢台'
      coordinate(41)%name='延安'
      coordinate(42)%name='平凉'
      coordinate(43)%name='锡林浩特'
      coordinate(44)%name='通辽'
      coordinate(45)%name='长春'
      coordinate(46)%name='赤峰'
      coordinate(47)%name='延吉'
      coordinate(48)%name='沈阳'
      coordinate(49)%name='临江'
      coordinate(50)%name='张家口'
      coordinate(51)%name='北京'
      coordinate(52)%name='乐亭'
      coordinate(53)%name='大连'
      coordinate(54)%name='章丘'
      coordinate(55)%name='荣成'
      coordinate(56)%name='青岛'
      coordinate(57)%name='那曲'
      coordinate(58)%name='拉萨'
      coordinate(59)%name='定和'
      coordinate(60)%name='托托河'
      coordinate(61)%name='玉树'
      coordinate(62)%name='达日'
      coordinate(63)%name='合作'
      coordinate(64)%name='武都'
      coordinate(65)%name='昌都'
      coordinate(66)%name='甘孜'
      coordinate(67)%name='红原'
      coordinate(68)%name='温江'
      coordinate(69)%name='巴塘'
      coordinate(70)%name='林芝'
      coordinate(71)%name='宜宾'
      coordinate(72)%name='西昌'
      coordinate(73)%name='丽江'
      coordinate(74)%name='咸宁'
      coordinate(75)%name='腾冲'
      coordinate(76)%name='昆明'
      coordinate(77)%name='思茅'
      coordinate(78)%name='蒙自'
      coordinate(79)%name='卢氏'
      coordinate(80)%name='郑州'
      coordinate(81)%name='汉中'
      coordinate(82)%name='泾河'
      coordinate(83)%name='南阳'
      coordinate(84)%name='安康'
      coordinate(85)%name='达县'
      coordinate(86)%name='恩施'
      coordinate(87)%name='宜昌'
      coordinate(88)%name='武汉'
      coordinate(89)%name='沙坪坝'
      coordinate(90)%name='长沙'
      coordinate(91)%name='怀化'
      coordinate(92)%name='贵阳'
      coordinate(93)%name='桂林'
      coordinate(94)%name='郴州'
      coordinate(95)%name='赣州'
      coordinate(96)%name='徐州'
      coordinate(97)%name='射阳'
      coordinate(98)%name='阜阳'
      coordinate(99)%name='南京'
      coordinate(100)%name='上海'
      coordinate(101)%name='安庆'
      coordinate(102)%name='杭州'
      coordinate(103)%name='南昌'
      coordinate(104)%name='衢州'
      coordinate(105)%name='洪家'
      coordinate(106)%name='邵武'
      coordinate(107)%name='福州'
      coordinate(108)%name='龙岩移动'
      coordinate(109)%name='河池'
      coordinate(110)%name='厦门'
      coordinate(111)%name='白色'
      coordinate(112)%name='梧州'
      coordinate(113)%name='清远'
      coordinate(114)%name='河源'
      coordinate(115)%name='汕头'
      coordinate(116)%name='南宁'
      coordinate(117)%name='北海'
      coordinate(118)%name='阳江'
      coordinate(119)%name='海口'
      coordinate(120)%name='三亚'
      coordinate(121)%name='西沙'


      coordinate(122)%stid='51848'; coordinate(122)%lambda=83.65; coordinate(122)%fi=37.93
      coordinate(123)%stid='52602'; coordinate(123)%lambda=93.38; coordinate(123)%fi=38.83
      coordinate(124)%stid='52889'; coordinate(124)%lambda=103.88; coordinate(124)%fi=36.05
      coordinate(125)%stid='54337'; coordinate(125)%lambda=121.12; coordinate(125)%fi=41.13
      coordinate(126)%stid='54374'; coordinate(126)%lambda=126.92; coordinate(126)%fi=41.80
      coordinate(127)%stid='54497'; coordinate(127)%lambda=124.33; coordinate(127)%fi=40.05
      coordinate(128)%stid='54776'; coordinate(128)%lambda=122.68; coordinate(128)%fi=37.40
      coordinate(129)%stid='54823'; coordinate(129)%lambda=116.98; coordinate(129)%fi=36.68
      coordinate(130)%stid='56294'; coordinate(130)%lambda=104.02; coordinate(130)%fi=30.67
      coordinate(131)%stid='56444'; coordinate(131)%lambda=98.90; coordinate(131)%fi=28.50
      coordinate(132)%stid='56989'; coordinate(132)%lambda=102.2; coordinate(132)%fi=23.01
      coordinate(133)%stid='57036'; coordinate(133)%lambda=108.93; coordinate(133)%fi=34.30
      coordinate(134)%stid='57290'; coordinate(134)%lambda=114.02; coordinate(134)%fi=33.00
      coordinate(135)%stid='57411'; coordinate(135)%lambda=106.80; coordinate(135)%fi=30.80
      coordinate(136)%stid='57515'; coordinate(136)%lambda=106.48; coordinate(136)%fi=29.52
      coordinate(137)%stid='58367'; coordinate(137)%lambda=121.43; coordinate(137)%fi=31.17
      coordinate(138)%stid='58666'; coordinate(138)%lambda=121.90; coordinate(138)%fi=28.45
      coordinate(139)%stid='58965'; coordinate(139)%lambda=121.50; coordinate(139)%fi=25.01
      coordinate(140)%stid='58968'; coordinate(140)%lambda=121.52; coordinate(140)%fi=25.03
      coordinate(141)%stid='59096'; coordinate(141)%lambda=114.40; coordinate(141)%fi=23.50
      coordinate(142)%stid='59287'; coordinate(142)%lambda=113.32; coordinate(142)%fi=23.13
      coordinate(143)%stid='59345'; coordinate(143)%lambda=119.57; coordinate(143)%fi=23.52
      coordinate(144)%stid='59362'; coordinate(144)%lambda=121.62; coordinate(144)%fi=24.02
      coordinate(145)%stid='59553'; coordinate(145)%lambda=120.43; coordinate(145)%fi=22.47
      coordinate(146)%stid='59647'; coordinate(146)%lambda=121.02; coordinate(146)%fi=22.30
      coordinate(147)%stid='59792'; coordinate(147)%lambda=116.72; coordinate(147)%fi=20.67
      coordinate(148)%stid='45004'; coordinate(148)%lambda=114.10; coordinate(148)%fi=22.50          !香港

      close(11)
end subroutine

subroutine char2num__3(a,num)

implicit none
      character*1 a(3)
      real num
      num=100*(ichar(a(1))-48)+10*(ichar(a(2))-48)+1*(ichar(a(3))-48)
end subroutine

subroutine char2num__2(a,num)
      implicit none
      character*1 a(2)
      real num
      num=10*(ichar(a(1))-48)+1*(ichar(a(2))-48)
end subroutine

function sixty2ten(hour,minute,second)
      implicit none
      real hour,minute,second
      real sixty2ten
      sixty2ten=hour+minute/60.0+second/3600.0
end function

end module

lqouc · 发表于 2014-4-30 20:02:20

这个程序是你自己写的？之前用过么？
我看主程序输出数据的地方好像是先输出经度的哈，应该先输出纬度，你换一下顺序试试看行不行。

风之叶 · 发表于 2014-4-30 20:32:23

lqouc 发表于 2014-4-30 20:02
这个程序是你自己写的？之前用过么？
我看主程序输出数据的地方好像是先输出经度的哈，应该先输出纬度，你 ...

谢谢啦!我去试试。这个程序是老师给的，我就直接拿来用了。

风之叶 · 发表于 2014-4-30 21:42:14

lqouc 发表于 2014-4-30 20:02
这个程序是你自己写的？之前用过么？
我看主程序输出数据的地方好像是先输出经度的哈，应该先输出纬度，你 ...

谢谢啦！非常感谢，我试了一下，图出来了，虽然不怎么美观，但也很不错了，谢谢学长的帮忙。

lqouc · 发表于 2014-4-30 21:48:12

风之叶发表于 2014-4-30 21:42
谢谢啦！非常感谢，我试了一下，图出来了，虽然不怎么美观，但也很不错了，谢谢学长的帮忙。

不用客气，欢迎常来论坛

步绝尘_南山居士 · 发表于 2014-5-1 00:38:16

这个一般是读入经纬度文件有问题造成的，你检查一下看看，有没有多读或者少读什么东西

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

[分享资料] DRADS插值时遇到warning from OACRES:less than two station.求指教！

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