上传一个自己编写的利用反距离权重插值函数

cnphoenix

        这些天正准备对缺测数据进行插值，经过清风的点拨，决定还是进行空间插值所得到的结果较好，于是编写了一个利用反距离权重法计算的点插值函数。
        欢迎讨论！

mofangbao

易小凯 发表于 2012-5-7 12:51

                               
登录/注册后可看大图

个人觉得用BARNES的插值比较方便一点。

有代码么
太好了，绝对支持原创~

陌小默的

支持，我童鞋也在 做，估计看到这个要乐了~

易小凯

个人觉得用BARNES的插值比较方便一点。

言深深

反距离权重是不是此前@topmad也做过？

topmad

言深深 发表于 2012-5-7 14:04

                               
登录/注册后可看大图

反距离权重是不是此前@topmad也做过？

恩，原来做过一次，再来学一下这位兄长的代码...

topmad

刚刚看了一下楼主的代码，应该是格点数据的插值方法吧  很清晰

cnphoenix

topmad 发表于 2012-5-7 16:30

                               
登录/注册后可看大图

刚刚看了一下楼主的代码，应该是格点数据的插值方法吧  很清晰

我准备用来对日值数据进行插值

易小凯

用BARNES方法将格点资料插值到站点资料上

1. DO it=1,NT
   
2. DO k=1,LEV
   
3.     DO index=1,64
   
4.         w=0.0
   
5.         DO i=1,NX1
   
6.         DO j=1,NY1
   
7.             r=(2.5*i+97.5-station(index).lon)**2+(2.5*j+12.5-station(index).lat)**2
   
8.             IF (r>6.25) THEN
   
9.                 CYCLE
   
10.             ELSE
    
11.                 station(index).ncep(k,it)=station(index).ncep(k,it)+exp(-r/6.25)*ncep(i,j,k,it)
    
12.                 w=w+exp(-r/6.25)
    
13.             ENDIF            
    
14.         ENDDO
    
15.         ENDDO
    
16.         station(index).ncep(k,it)=station(index).ncep(k,it)/w
    
17.     ENDDO
    
18. ENDDO
    
19. ENDDO
    
20. WRITE(*,*) 'NCEP INTERPOLATION OK!'
    
21. 
    

DO it=1,NT<br /> DO k=1,LEV<br />     DO index=1,64<br />         w=0.0<br />         DO i=1,NX1<br />         DO j=1,NY1<br />             r=(2.5*i+97.5-station(index).lon)**2+(2.5*j+12.5-station(index).lat)**2<br />             IF (r>6.25) THEN<br />                 CYCLE<br />             ELSE<br />                 station(index).ncep(k,it)=station(index).ncep(k,it)+exp(-r/6.25)*ncep(i,j,k,it)<br />                 w=w+exp(-r/6.25)<br />             ENDIF            <br />         ENDDO<br />         ENDDO<br />         station(index).ncep(k,it)=station(index).ncep(k,it)/w<br />     ENDDO<br /> ENDDO<br /> ENDDO<br /> WRITE(*,*) 'NCEP INTERPOLATION OK!'<br /> <br />

因为格点分辨率是2.5*2.5，故选6.25为阈值，对应的数据分辨率改一下阈值就好。
阈值大了，会滤去一些信息；阈值小了，可能有效半径内没有格点值

cnphoenix

易小凯 发表于 2012-5-7 18:48

                               
登录/注册后可看大图

用BARNES方法将格点资料插值到站点资料上

对于栅格点上的插值不懂，学习了