降水图画出来跟底图不匹配

Mr.Yes

如图1，降水图画出来跟底图不匹配，我是用micaps里的第三类数据画的， 原始资料如图2，下面的才是降水吧，从左到右分别是站号，经度，纬度，海拔，降水量。不是很清楚上面一坨数据是什么意思。
我先将站点资料处理为二进制文件：
以下是fortran程序：

  

1. !将micaps数据格式第3类降水数据转化为GRADS可读格式
   
2.   parameter(nzd=2000)
   
3.   character*8 :: stid(nzd)
   
4.   integer :: k,zd
   
5.   real :: r(nzd),rlon(nzd),rlat(nzd),hgt(nzd)
   
6.   character*100 :: chf1,chf2,chf3,chf4,chf5,chf6,chf7,chf8,chf9,chf10,chf11,chf12,chf13,chf14,chf15
   
7.   character*30 :: cfile='G:\data\r24-8-p\'
   
8.   character*30 :: cfile1='14061708.000'
   
9. !-----------------------------------------
   
10.        open(1,file=trim(cfile)//trim(cfile1),status='old')
    
11.        read(1,500)chf1
    
12.        read(1,500)chf2
    
13.        read(1,500)chf3
    
14.        read(1,500)chf4
    
15.     read(1,500)chf5
    
16.        read(1,500)chf6
    
17.        read(1,500)chf7
    
18.        read(1,500)chf8
    
19.     read(1,500)chf9
    
20.        read(1,500)chf10
    
21.        read(1,500)chf11
    
22.     read(1,500)chf13
    
23.     read(1,500)chf14
    
24.        read(1,500)chf15
    
25. 
    
26. 500 format(a100)
    
27.      !以上是把站号以上的数据读掉
    
28.     read(chf15,*)k,zd  !读这个文件里面有多少个站，一共有zd个站，k是把前面的1读掉
    
29.      do i=1,zd
    
30.         read(1,*)stid(i),RLon(i),RLat(i),hgt(i),R(i)
    
31.      enddo
    
32. !以上程序已经把站号，经纬度，海拔，降水量都读取完毕，下面把数据写到2进制文件中
    
33.      open(2,file='G:\data\out\6\rain.grd',form='binary')
    
34.   tim=0.0
    
35.   nlev=1
    
36.   nflag=1
    
37.   do 40 i=1,zd
    
38.   WRITE(2)stid(i),rlat(i),rlon(i),tim,nlev,nflag,r(i)
    
39.      40 continue
    
40.   nlev=0
    
41.   WRITE(2)stid(i-1),rlat(i-1),rlon(i-1),tim,nlev,nflag
    
42.   CLOSE(2)
    
43.      print*,stid(5),RLon(5),RLat(5),hgt(5),R(5)
    
44.   end
    
45. 
    
46. !***********************************************      
    
47. !此程序对应rain.ctl文件，运行！ stnmap+rain.ctl生成*.map文件
    

!将micaps数据格式第3类降水数据转化为GRADS可读格式<br />   parameter(nzd=2000)<br />   character*8 :: stid(nzd)<br />   integer :: k,zd<br />   real :: r(nzd),rlon(nzd),rlat(nzd),hgt(nzd)<br />   character*100 :: chf1,chf2,chf3,chf4,chf5,chf6,chf7,chf8,chf9,chf10,chf11,chf12,chf13,chf14,chf15<br />   character*30 :: cfile='G:\data\r24-8-p\'<br />   character*30 :: cfile1='14061708.000'<br /> !-----------------------------------------<br />        open(1,file=trim(cfile)//trim(cfile1),status='old')<br />        read(1,500)chf1<br />        read(1,500)chf2<br />        read(1,500)chf3<br />        read(1,500)chf4<br />     read(1,500)chf5<br />        read(1,500)chf6<br />        read(1,500)chf7<br />        read(1,500)chf8<br />     read(1,500)chf9<br />        read(1,500)chf10<br />        read(1,500)chf11<br />     read(1,500)chf13<br />     read(1,500)chf14<br />        read(1,500)chf15<br /> <br /> 500 format(a100)<br />      !以上是把站号以上的数据读掉<br />     read(chf15,*)k,zd  !读这个文件里面有多少个站，一共有zd个站，k是把前面的1读掉<br />      do i=1,zd<br />         read(1,*)stid(i),RLon(i),RLat(i),hgt(i),R(i)<br />      enddo<br /> !以上程序已经把站号，经纬度，海拔，降水量都读取完毕，下面把数据写到2进制文件中<br />      open(2,file='G:\data\out\6\rain.grd',form='binary')<br />   tim=0.0<br />   nlev=1<br />   nflag=1<br />   do 40 i=1,zd<br />   WRITE(2)stid(i),rlat(i),rlon(i),tim,nlev,nflag,r(i)<br />      40 continue<br />   nlev=0<br />   WRITE(2)stid(i-1),rlat(i-1),rlon(i-1),tim,nlev,nflag<br />   CLOSE(2)<br />      print*,stid(5),RLon(5),RLat(5),hgt(5),R(5)<br />   end<br /> <br /> !***********************************************      <br /> !此程序对应rain.ctl文件，运行！ stnmap+rain.ctl生成*.map文件<br />

与rain.grd相应的ctl文件:
dset g:\data\out\6\rain.grd
dtype station
stnmap g:\data\out\rain1.map
undef -999.0
title rain
tdef 1 linear 08:00Z17JUL2014 24hr
vars 1
p 0 99 rain
endvars

然后生成格点文件，以下是fortran程序：

1. !生成格点文件
   
2. !grid.grd文件的经度要高于或等于rain.grd的精度（取1度就好）
   
3. !grid.grd文件的范围要大于或等于rain.grd的范围。
   
4. !grid.grd文件的每个点上均赋值为1
   
5. !grid.grd文件的描述文件中时间说明一定要与rain.ctl中时间一致。
   
6. parameter(nx=71,ny=51)
   
7. real lat(ny),lon(nx)
   
8. real s(nx,ny)
   
9. open(1,file='G:\data\out\6\grid.grd',form='binary')
   
10. lat(1)=5.0
    
11. lon(1)=70
    
12. do j=1,ny-1
    
13. lat(j+1)=lat(j)+1.0
    
14. enddo
    
15. do i=1,nx-1
    
16. lon(i+1)=lon(i)+1.0
    
17. enddo
    
18.   do i=1,nx
    
19.   do j=1,ny
    
20.   s(i,j)=1
    
21.   enddo
    
22.   enddo
    
23. write(1)s
    
24. end
    

复制代码

与格点数据grid.grd对应的ctl文件：
dset G:\data\out\6\grid.grd
undef -999.0
title Sample GRID Data
xdef 71 linear 70 1
ydef 41 linear 15 1
zdef 1  linear 500 1
tdef 1  linear 08:00Z17JUL2014 24hr
vars 1
g 0 99 grid data
endvars

最后是画图的gs文件：
'reinit'
'open G:\data\out\6\grid.ctl'
'open G:\data\out\6\rain.ctl'
'set grads off'
'enable print G:\data\out\6\picture\rain.gmf'
'set lon 70 135'
'set lat 10 55'
'set mpdset cnworld'
'define a=oacres(g,p.2,10,7,4,2,1)'
'define a1=maskout(a,g-0.5)'
'define aa=smth9(a1)'
'set xlopts 1 10 0.18'
'set ylopts 1 10 0.18'
'd aa'
'print'
'disable print'
;
这里用到了cnworld的底图，网上都有下载的，结果出来就成了之前的样子，不知道什么地方出了问题，恳请各位前辈予以解答，谢谢！

river

站点资料，利用插值才能画出等值线。既然是插值，就不可能正正好好插值到你要的边界上，有超出国界范围的才是正常的。你需要做的就是把国界外的覆盖掉就可以了

yuenai

如果经纬度没对应错的话  用个mask文件把边缘盖掉应该就好了

545881287

我也出现过同类问题。