网格点插值站点循环多个时次

动动

此程序为网格点数据插值到站点，一个时次。
请大神指点在哪里加循环，改为4个时次的循环

1. 
   
2. !   双线性插值，格点到站点，考虑权重  !
   
3. !      modifioed from    mofangbao(清风) by djj!   
   
4. 
   
5. !===================================================!
   
6. program bilineargird2stn
   
7. integer,parameter::oldx=97,oldy=61    !格点数据xy向网格点数
   
8. integer,parameter::newxy=1308          !站点数
   
9. integer,parameter::time=4    !时间层次
   
10. real,parameter::oldLon_start=100.0,oldLat_start=20.0,oldInterpx=0.25,oldInterpy=0.25  !格点数据起始经纬度和格距
    
11. real oldVal(oldx,oldy),newVal(newxy),newPoint(newxy,2)  !分别是：格点数据，站点数据，站点数据经纬度，newVal的点顺序和nowPoint中的对应
    
12. character*8 ::stid
    
13. integer::i,j,k
    
14. logical isFirst   !判断是否第一个时次
    
15. isFirst=.true.      !后面无需手动修改此值，程序会自动修改
    
16. 
    
17. 
    
18. !!读取站点经纬度
    
19. open(27,file='C:\PCGrADS\win32\st1308-2.txt')
    
20. do i=1,newxy
    
21.     read(27,*) stid,newPoint(i,2),newPoint(i,1)  !先经度后纬度
    
22. enddo
    
23. close(27)
    
24. 
    
25. 
    
26. !!读取网格点数据
    
27. open(29,file='C:\PCGrADS\win32\24cma2015081619.grib',form='binary')
    
28. read(29)  ((oldVal(i,j),i=1,oldx),j=1,oldy)
    
29. close(29)
    
30. 
    
31. 
    
32. !do k=1,4
    
33. !!调用插值程序,这里只哟一个时次，因此没有循环
    
34. call getNewVal(oldVal,newVal,oldLon_start,oldLat_start,oldInterpx,oldInterpy,newPoint,isFirst)
    
35. open(25,file='C:\PCGrADS\win32\grid2stn.txt')
    
36. do i=1,newxy
    
37. write(25,*) newVal(i)
    
38. enddo
    
39. end
    
40. 
    
41. 
    
42. 
    
43. 
    
44. subroutine getNewVal(oldVal,newVal,oldLon_start,oldLat_start,oldInterpx,oldInterpy,newPoint,isFirst)
    
45. integer,parameter::oldx=97,oldy=61,newxy=1308
    
46. real::oldLon_strat,oldLat_strat
    
47. real,parameter::undef=-999.0  
    
48. real::oldInterpx,oldInterpy      
    
49. integer pos(newxy,4),nstart,nend  !weight 1,2是左侧点和右侧点的分别权重，3,4是下面点和上面点分别的权重
    
50. real tmpVal_up,tmpVal_down !本次纬向下方和纬向上方的线性插值结果以及权重
    
51. real::weight(newxy,4)
    
52. real newLon,newLat,weight_end  !本次插值的新点位置，以及返回的权重
    
53. real oldVal(oldx,oldy),newVal(newxy),newPoint(newxy,2)  !旧值，新值，新点经纬度
    
54. logical isfirst  !是否第一个时次，只有第一个时次需要计算四周格点分布以及权重，后面的套用即可
    
55. 
    
56. if(isfirst)then
    
57. !获取新点的四周的点以及每个点的权重
    
58. do i=1,newxy
    
59.   newLon=newPoint(i,1)
    
60. 
    
61. ! print*,newLon,oldLon_start
    
62.   call getPosAndWeight(newLon,oldLon_start,nstart,nend,weight_end,oldInterpx,isFirst)
    
63.   !print*,newLon,nstart,nend,weight_end
    
64.   if(nend>oldx)then
    
65.     if(nstart==oldx)then
    
66.       pos(i,1)=nstart
    
67.       pos(i,2)=nstart
    
68.       weight(i,1)=1
    
69.       weight(i,2)=0
    
70.     else
    
71.       pos(i,1)=-1
    
72.       pos(i,2)=-1  !超出了原来格点的范围，标记为无效值
    
73.       weight(i,1)=-1
    
74.       weight(i,2)=-1
    
75.     endif
    
76.   else
    
77.     pos(i,1)=nstart
    
78.     pos(i,2)=nend
    
79.     weight(i,1)=1-weight_end
    
80.     weight(i,2)=weight_end
    
81.   endif
    
82. 
    
83. !print*, pos(i,1),pos(i,2)
    
84. !write(*,*) weight(i,1),weight(i,2)
    
85. 
    
86. 
    
87.   newLat=newPoint(i,2)
    
88. ! write(*,*) newLat,oldLat_start
    
89.   call getPosAndWeight(newLat,oldLat_start,nstart,nend,weight_end,oldInterpy,isFirst)
    
90. 
    
91.   if(nend>oldy)then
    
92.     if(nstart==oldy)then
    
93.       pos(i,3)=nstart
    
94.       pos(i,4)=nstart
    
95.       weight(i,3)=1
    
96.       weight(i,4)=0
    
97.     else
    
98.       pos(i,3)=-1
    
99.       pos(i,4)=-1  !超出了原来格点的范围，标记为无效值
    
100.       weight(i,3)=-1
     
101.       weight(i,4)=-1
     
102.     endif
     
103.   else
     
104.     pos(i,3)=nstart
     
105.     pos(i,4)=nend
     
106.     weight(i,3)=1-weight_end
     
107.     weight(i,4)=weight_end
     
108.   endif
     
109. 
     
110. !print*, pos(i,3),pos(i,4)
     
111. !write(*,*) weight(i,3),weight(i,4)  
     
112. !pause   
     
113. enddo
     
114. endif
     
115. 
     
116. 
     
117. 
     
118. !开始插值
     
119. do i=1,1308
     
120.   if(pos(i,1)<0 .or.pos(i,3)<0)then
     
121.     newval(i)=-999.0
     
122.   else
     
123.     !先做纬向插值
     
124.     if(pos(i,1)==undef .or. pos(i,2)==undef .or. pos(i,3)==undef .or. pos(i,4)==undef)then
     
125.       newVal(i)=-9999.0
     
126.     else
     
127.       if(oldVal(pos(i,1),pos(i,3))>0.and.oldVal(pos(i,2),pos(i,3))>0) tmpVal_down=weight(i,1)*oldVal(pos(i,1),pos(i,3))+weight(i,2)*oldVal(pos(i,2),pos(i,3))
     
128.       if(oldVal(pos(i,1),pos(i,3))>0.and.oldVal(pos(i,2),pos(i,3))<0) tmpVal_down=oldVal(pos(i,1),pos(i,3))
     
129.       if(oldVal(pos(i,1),pos(i,3))<0.and.oldVal(pos(i,2),pos(i,3))>0) tmpVal_down=oldVal(pos(i,2),pos(i,3))
     
130.       if(oldVal(pos(i,1),pos(i,3))<0.and.oldVal(pos(i,2),pos(i,3))<0) tmpVal_down=-999.0
     
131.       if(oldVal(pos(i,1),pos(i,4))>0.and.oldVal(pos(i,2),pos(i,4))>0) tmpVal_up=weight(i,1)*oldVal(pos(i,1),pos(i,4))+weight(i,2)*oldVal(pos(i,2),pos(i,4))
     
132.       if(oldVal(pos(i,1),pos(i,4))<0.and.oldVal(pos(i,2),pos(i,4))>0) tmpVal_up=oldVal(pos(i,2),pos(i,4))
     
133.       if(oldVal(pos(i,1),pos(i,4))>0.and.oldVal(pos(i,2),pos(i,4))<0) tmpVal_up=oldVal(pos(i,1),pos(i,4))
     
134.       if(oldVal(pos(i,1),pos(i,4))<0.and.oldVal(pos(i,2),pos(i,4))<0) tmpVal_up=-999.0
     
135.       if(tmpVal_down>0.and.tmpVal_up>0) then
     
136.         newVal(i)=weight(i,3)*tmpVal_down+weight(i,4)*tmpVal_up
     
137.       else if(tmpVal_down>0.and.tmpVal_up<0) then
     
138.         newVal(i)=tmpVal_down
     
139.       else if(tmpVal_down<0.and.tmpVal_up>0) then
     
140.         newVal(i)=tmpVal_up
     
141.       else
     
142.         newVal(i)=-999.0
     
143.       endif
     
144.     endif
     
145.   endif
     
146. enddo
     
147. 
     
148. 
     
149. open(25,file='C:\PCGrADS\win32\weight.txt')
     
150. do i=1,1308
     
151. write(25,*) weight(i,:)
     
152. enddo
     
153. close(25)
     
154. endsubroutine
     
155. 
     
156. !计算某个方向的两侧的点以及后侧点（上方点）的权重
     
157. subroutine  getPosAndWeight(newPos,oldstart,nstart,nend,weight,oldInterp,isfirst)
     
158. integer nstart,nend  !两侧点
     
159. real newPos,weight,oldInterp
     
160. real oldstart   
     
161. logical isfirst
     
162. !离的近的点权重大，但是占的距离比重小，所以权重=1-距离比重，经纬度从0开始，但是标号从1开始，所以加1
     
163. print*,newPos,oldstart,oldInterp
     
164. nstart=int((newPos-oldstart)/oldInterp)+1
     
165. nend=nstart+1
     
166. weight=(newPos-(nstart-1)*oldInterp-oldstart)/oldInterp  !  |---|-------|  此表达式计算的是后一段的权重（前一段的距离比重），该段总长为interp
     
167. isfirst=.false.
     
168. endsubroutine

复制代码

mofangbao

这个。。。让别人花钱给你解决问题，也是醉了。。。如果楼主是想奖励帮助你的人，应该发悬赏帖

动动

mofangbao 发表于 2016-1-11 20:57
这个。。。让别人花钱给你解决问题，也是醉了。。。如果楼主是想奖励帮助你的人，应该发悬赏帖

原谅我吧，第一次发帖。我是想奖励帮助我的人，可是也不知道怎么弄。还麻烦大神帮我看一下，谢谢

zhangy-sjyy

没数据，没办法进行

小娜子

不错   支持！

弄晴

楼主你好，最近也在做格点插值站点，请问问题解决了吗

aridlhy

加一个时间循环吧