按照清风大大的grads站点文件.pdf中的站点插值到格点，最后出图有误

xmgcqlmw

1

现有1979年6月的气象台753站资料，想把站点插值到格点，画1979年6月1日降水图，站点数据按经纬度输出画的降水量分布图没有问题，但是最终插值后画的等值线图却不对，不知道是哪里出了问题，请问有谁知道吗？？感激不尽
  以下是我的txt以及所用到的程序：
我的txt格式如下：分别是站点号，纬度，经度，年，月，日，降水量0.1mm

1. 54423 40.59 117.57 1979 6 13 0.4
   
2. 54423 40.59 117.57 1979 6 14 16.0
   
3. 54423 40.59 117.57 1979 6 15 0.1
   
4. 54423 40.59 117.57 1979 6 16 0.4
   
5. 54423 40.59 117.57 1979 6 17 0.0
   
6. 54423 40.59 117.57 1979 6 18 0.0
   
7. 54423 40.59 117.57 1979 6 19 8.9
   
8. 54423 40.59 117.57 1979 6 20 1.0
   
9. 54423 40.59 117.57 1979 6 21 99999.0
   
10. 54423 40.59 117.57 1979 6 22 0.0
    
11. 54423 40.59 117.57 1979 6 23 99999.0
    
12. 54423 40.59 117.57 1979 6 24 36.3
    
13. 54423 40.59 117.57 1979 6 25 21.7
    
14. 54423 40.59 117.57 1979 6 26 99999.0
    
15. 54423 40.59 117.57 1979 6 27 0.0
    
16. 54423 40.59 117.57 1979 6 28 0.0
    
17. 54423 40.59 117.57 1979 6 29 3.1
    
18. 54423 40.59 117.57 1979 6 30 4.5
    
19. 54429 40.12 117.57 1979 6 1 0.0
    
20. 54429 40.12 117.57 1979 6 2 4.5
    
21. 54429 40.12 117.57 1979 6 3 99999.0
    
22. 54429 40.12 117.57 1979 6 4 0.0
    
23. 54429 40.12 117.57 1979 6 5 60.5
    
24. 54429 40.12 117.57 1979 6 6 7.5
    
25. 54429 40.12 117.57 1979 6 7 11.7

复制代码

2、 在 fortran 中读入这些资料并按照 GrADS 的规定将数据进行输出

1. Program sta2grd
   
2. Implicit none
   
3. !这里是程序的变量声明
   
4. Character*8 stid
   
5. Real lon,lat,year,month,day,rain,tim
   
6. integer nlev,flag
   
7. !变量声明结束
   
8. !程序开始
   
9. tim=0.0
   
10. nlev=1
    
11. flag=1
    
12. Open(1,file='F:\jiangshui\753\1979-2014\out\oldrain\pipei\1979_6pipei.txt',status='old')
    
13. open(2,file='F:\jiangshui\753\1979-2014\out\oldrain\grd\1979601sta.grd',status='replace',form='binary')
    
14. 10 Read(1,*,end=100)stid,lat,lon,year,month,day,rain
    
15. if(day==1.0)then
    
16. 
    
17. !缺测处理
    
18. if(int(rain)==99999)then   
    
19. rain=99999
    
20. else
    
21. rain=rain*10
    
22. endif
    
23. Print*,stid,lat,lon,rain,day
    
24. write(2)stid,lat,lon,tim,nlev,flag,rain
    
25. 
    
26. else
    
27. end if
    
28. !pause
    
29. Goto 10
    
30. 100 continue
    
31. Close(1)
    
32. nlev=0
    
33. 
    
34. Write(2)stid,lat,lon,tim,nlev,flag
    
35. 
    
36. close(2)
    
37. !程序结束
    
38. End

复制代码

3、 为生成的站点数据编写 CTL（数据描述文件）

1. DSET F:/jiangshui/753/1979-2014/out/oldrain/grd/1979601sta.grd
   
2. DTYPE station
   
3. STNMAP F:/jiangshui/753/1979-2014/out/oldrain/grd/sta.map
   
4. UNDEF 99999.00
   
5. TITLE Rain Data Sample
   
6. *这里的时间可以根据上面示例文件中的第一行的时间填写
   
7. TDEF 1 linear 01jun1979 24hr
   
8. VARS 1
   
9. rain 0 99 Rain Data
   
10. ENDVARS

复制代码

4、 为站点的二进制数据生成站点映射文件（ *.map 文件）      有生成成功
5. 看看站点数据如何
     画的图如图1

6、 生成一个格点的背景场，为插值做准备

1. Program gen_grid
   
2. Implicit none
   
3. !这里是定义的变量
   
4. integer i,j
   
5. integer,parameter::x=29,y=17
   
6. real grid
   
7. !程序开始
   
8. grid=1.0
   
9. Open(1,file='F:\jiangshui\753\1979-2014\out\oldrain\grd\1979601grid.grd',status='replace',form='binary')
   
10. Do i=1,y
    
11. Do j=1,x
    
12. Write(1)grid
    
13. Enddo
    
14. Enddo
    
15. close(1)
    
16. !程序结束
    
17. End

复制代码

7、为生成的格点场配置 CTL 文件

1. DSET F:\jiangshui\753\1979-2014\out\oldrain\grd\1979601grid.grd
   
2. TITLE Grid Data Sample
   
3. UNDEF -999.0
   
4. XDEF 29 LINEAR 70.0 2.5
   
5. YDEF 17 LINEAR 15.0 2.5
   
6. ZDEF 1 LEVELS 1000
   
7. TDEF 1 LINEAR 01jul1979 24hr
   
8. VARS 1
   
9. g 0 99 Grid Data
   
10. ENDVARS

复制代码

8、编写 gs 脚本进行插值后作出图形

1. 'reinit'
   
2. 'open F:\jiangshui\753\1979-2014\out\oldrain\grd\grid.ctl'
   
3. 'open F:\jiangshui\753\1979-2014\out\oldrain\grd\sta.ctl'
   
4. 'set lon 70 140'
   
5. 'set lat 15 55'
   
6. 'set mpdset hires'
   
7. 
   
8. 'define a=oacres(g,rain.2,8,5,3,2,1)'
   
9. 'define al=maskout(a,g-1)'
   
10. 'define aa=smth9(al)'
    
11. 
    
12. 'd aa'
    
13. 'printim F:\jiangshui\753\1979-2014\out\oldrain\grd\sta_grid32.png white'
    
14. ;
    

复制代码

最后出的图（图2）值特别大，跟站点分布图的值不一样，研究了很久，实在不知道哪里错了，希望各位大大提点，不甚感激。

river

原始资料里的降水的缺测值全部都处理成99999.0了？我看第一张图上还有30开头的五位数，那个应该没有处理吧

xmgcqlmw

附上清风的教程，向大大们学习

男紫汗

没考虑缺测值

xmgcqlmw

男紫汗 发表于 2017-7-4 19:45
没考虑缺测值

我缺测值设置了99999，也试过了-9.99E33，但是还是不对

男紫汗

第一张图看着不正常吧，日降水能有上千？

xmgcqlmw

river 发表于 2017-7-5 12:06
原始资料里的降水的缺测值全部都处理成99999.0了？我看第一张图上还有30开头的五位数，那个应该没有处理吧

确实是有，之前都没注意到，谢谢，我重新排查一遍

xmgcqlmw

男紫汗 发表于 2017-7-4 20:52
第一张图看着不正常吧，日降水能有上千？

是的是的，有些数据好像没处理成功，我重新排查一遍，谢谢

xmgcqlmw

男紫汗 发表于 2017-7-4 20:52
第一张图看着不正常吧，日降水能有上千？

原始资料里1979年6月24日降水量（单位0.1mm）是1472，请问日降水量147.2mm这个值可信吗，能用还是要做处理

男紫汗

xmgcqlmw 发表于 2017-7-5 13:09
原始资料里1979年6月24日降水量（单位0.1mm）是1472，请问日降水量147.2mm这个值可信吗，能用还是要做处 ...

能用。。。。