[已解決]有關站點資料畫等值線圖問題想請教

addname · 发表于 2016-3-21 14:46:50

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

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

x

本帖最后由 addname 于 2016-3-24 14:01 编辑

想請教一下，讀取完站點資料後，將資料點畫圖出來很正常，

站點資料

格點資料也可以正常畫圖，但是插值畫等值線圖時卻變成底下這模樣，不知道是哪裡出了問題

等值線圖

station.txt (13.5 KB, 下载次数: 102)

2016-3-21 14:36 上传

点击文件名下载附件
經緯度

201511_auto_hr.txt (14.71 MB, 下载次数: 296)

2016-3-21 14:35 上传

点击文件名下载附件
站點資料
----------------------------程式碼--------------------------------------------------
Text1.f90
program test
Implicit none
Character*8 stid,id
Real lon,lat,time,temp
integer nlev,nflag,i,j,ii,jj,it
integer date,dd,hh

open(21,file='201511.dat',status='unknown',FORM='UNFORMATTED',ACCESS='STREAM')
   do it=1,720
   dd=int(it/24)
      hh=it-dd*24

      if(hh .eq. 0)then
      dd=int(it/24)-1
      hh=24
end if
   do i=1,512
open(11, file='station.txt', form='formatted', status='old')
read(11,*,end=90)stid,lon,lat
time=0.
nlev=1
nflag=1
open(12,file='201511_auto_hr.txt',form='formatted',status='old')
   do j=1,226800
read(12,'(a6,1x,i10,10x,f4.1)',end=92) id,date,temp
   if(id .eq. stid .and. date .eq. (2015110100+dd*100+hh)) then
write(21) stid,lat,lon,time,nlev,nflag
write(21) temp
   go to 91
end if
   end do
92 close(12)
write(21) stid,lat,lon,time,nlev,nflag
write(21) 99.99
91 end do
90 close(11)
nlev=0
write(21)stid,lat,lon,time,nlev,nflag
end do
stop
end

st.ctl
dset ^201511.dat
dtype  station
stnmap ^201511.map
undef  99.99
title  Station Data Sample
tdef 720 linear 01z01nov2015 1hr
vars 1
temp  0  99  temp
endvars

st.gs
'reinit'
'open st.ctl'
'set grads off'
'set lon 119.8 122'
'set lat 21.9 25.3'
'set t 8'
'set mpdset hires'
'set dignum 1'
'd temp'
'draw title 2015110108 temp'
'printim 110108_temp.png white png'
'quit'

Text2.f90
program test2

Implicit none

integer i,j,t
real lat(440),lon(220),grid(220,440,720)


open(1,file='grid.grd',form='unformatted',access='direct',recl=220*440*4)
lat(1)=21.9
lon(1)=119.8

do j=1,440-1
lat(j+1)=lat(j)+0.01
end do

do i=1,220-1
lon(i+1)=lon(i)+0.01
end do

grid=99.99

do t=1,720
      write(1,rec=t)((grid(i,j,t),i=1,220),j=1,440)
      end do

end

test.ctl
dset grid.grd
title  Station Data Sample
undef 99.99
xdef 220 linear 119.8 0.01
ydef 440 linear 21.9  0.01
zdef 1 levels 1000
tdef  720 linear 01z01nov2015 1hr
vars 1
g  0  99  Grid Data
endvars

test.gs
'reinit'
'open test.ctl'
'open st.ctl'
'set lon 119.8 122'
'set lat 21.9 25.3'
'set t 8'
'set mpdset hires'
'set gxout shaded'
'set gxout contour'
'set dignum 1'
'd oacres(g,temp.2)'
'printim test.png white png'
'quit'

MeteoInfo · 发表于 2016-3-22 15:55:18

addname 发表于 2016-3-22 09:36
謝謝您，已經可以畫出來了，原來是我的分辨率太大所以畫不出來，一直以為格點數要多，才夠多格 ...

凑个热闹：MeteoInfoLab脚本示例：站点数据读取插值http://bbs.06climate.com/forum.p ... 091&fromuid=106
(出处: 气象家园)

兰溪之水 · 发表于 2016-3-21 18:40:29

测试了下，你生成二进制文件的效率太低了，在你的基础上，重新改了下，还有你原来的程序读取温度的数据是有问题的，因为你的温度数据存在-9997的值，猜测可能是缺测值，所以你这样读“read(12,'(a6,1x,i10,10x,f4.1)') id,date,temp”不是一个robust的做法，给你做了相应的修改~具体看

program test
Implicit none
integer,parameter :: nstid = 512, nrecord = 226800
Character*8 :: stid(nstid), id(nrecord)
Real :: lon(nstid),lat(nstid),time,slp,temp(nrecord)
integer :: nlev,nflag,i,j,it
integer :: date(nrecord),dd,hh
open(11, file='station.txt', form='formatted', status='old')
do i=1, nstid
read(11,*)stid(i),lon(i),lat(i)
!write(*,*) stid(i),lon(i),lat(i)
end do
open(12,file='201511_auto_hr.txt',form='formatted',status='old')
do j=1, nrecord
!read(12,'(a6,1x,i10,10x,f4.1)') id(j),date(j),temp(j)
read(12,*) id(j),date(j),slp,temp(j)
if(abs(temp(j)+9997).lt.1.e-6) temp(j) = 99.99
end do
open(21,file='201511.dat',status='unknown',FORM='UNFORMATTED',ACCESS='STREAM')
do it=1,720
write(*,*) "time: ",it
dd=int(it/24)
hh=it-dd*24
if(hh .eq. 0)then
dd=int(it/24)-1
hh=24
end if
time=0.
nlev=1
nflag=1
do i=1, nstid
do j=1, nrecord
if(id(j) .eq. stid(i) .and. date(j) .eq. (2015110100+dd*100+hh)) then
write(21) id(j),lat(i),lon(i),time,nlev,nflag,temp(j)
exit
end if
end do
end do
nlev=0
write(21) stid(nstid),lat(nstid),lon(nstid),time,nlev,nflag
end do
close(21)
stop
end

program test 
    Implicit none 
    integer,parameter :: nstid = 512, nrecord = 226800 
    Character*8       :: stid(nstid), id(nrecord) 
    Real              :: lon(nstid),lat(nstid),time,slp,temp(nrecord) 
    integer           :: nlev,nflag,i,j,it 
    integer           :: date(nrecord),dd,hh 
    
    open(11, file='station.txt', form='formatted', status='old') 
    do i=1, nstid 
        read(11,*)stid(i),lon(i),lat(i) 
        !write(*,*) stid(i),lon(i),lat(i) 
    end do 
    open(12,file='201511_auto_hr.txt',form='formatted',status='old') 
    do j=1, nrecord 
        !read(12,'(a6,1x,i10,10x,f4.1)') id(j),date(j),temp(j) 
        read(12,*) id(j),date(j),slp,temp(j) 
        if(abs(temp(j)+9997).lt.1.e-6) temp(j) = 99.99 
    end do 
      
    open(21,file='201511.dat',status='unknown',FORM='UNFORMATTED',ACCESS='STREAM') 
    do it=1,720 
        write(*,*) "time: ",it 
         
        dd=int(it/24) 
        hh=it-dd*24 
        if(hh .eq. 0)then 
            dd=int(it/24)-1 
            hh=24 
        end if 
    
        time=0. 
        nlev=1 
        nflag=1 
        do i=1, nstid 
            do j=1, nrecord 
                if(id(j) .eq. stid(i) .and. date(j) .eq. (2015110100+dd*100+hh)) then 
                    write(21) id(j),lat(i),lon(i),time,nlev,nflag,temp(j) 
                    exit 
                end if 
            end do 
        end do 
        nlev=0 
        write(21) stid(nstid),lat(nstid),lon(nstid),time,nlev,nflag 
    end do 
    close(21) 
    
    stop 
    end

插值成网格数据是因为你的背景的网格分辨率太高（也就是格距太小），然后插值半径默认又太小，所以才会出现这样的结果~
可以把插值半径改大，比如'd oacres(g,temp.2,50,40,30,20,10,5,1)'，当然我这是夸大了，或者你把背景的网格分辨率变粗一点，
比如把0.01改成0.1，即
xdef 220 linear 119.8 0.1
ydef 440 linear 21.9 0.1
当然最后还是自己慢慢调，自己觉得满意即可~Good luck！
@传说中的谁 @mofangbao 这活可好

传说中的谁 · 发表于 2016-3-21 15:16:37

可能的原因：一是插值背景网格的分辨率，二是插值的影响半径，即oacres函数后面的参数@兰溪之水

addname · 发表于 2016-3-21 17:33:27

嗯，應該不是oacres函数后面的参数的問題，之前有修改過，結果也都差不多。
所以是分辨率嚕？您的意思是格點數的問題嗎？我再將每個格點的間隔變小看看，謝謝！

weiqiyun0049 · 发表于 2016-3-21 17:55:56

本帖最后由 weiqiyun0049 于 2016-3-21 17:57 编辑

addname 发表于 2016-3-21 17:33
嗯，應該不是oacres函数后面的参数的問題，之前有修改過，結果也都差不多。
所以是分辨率嚕？您的意思是格 ...

不是变小吧。。。变大会好点，或者你用gxout shaded看看是什么鬼就好了，把gxout contour注释掉

兰溪之水 · 发表于 2016-3-21 18:43:57

简单测试图

syl136 · 发表于 2016-3-21 18:59:20

我是来看评论的

李璠 · 发表于 2016-3-22 09:12:14

看到大神出现很开心

addname · 发表于 2016-3-22 09:36:21

本帖最后由 addname 于 2016-3-22 10:02 编辑

謝謝您，已經可以畫出來了，原來是我的分辨率太大所以畫不出來，一直以為格點數要多，才夠多格點來畫等值線，結果是相反的。
真是太感激您了，如果有問題再請教您。

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

[秀图] [已解決]有關站點資料畫等值線圖問題想請教

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

相关帖子