分享我写FORTRAN程序（反距离权重法插值）的一些步骤和心得

咕噜 · 发表于 2016-1-16 18:51:17

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

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

x

本帖最后由咕噜于 2018-7-1 11:40 编辑

想要把再分析格点资料插值到站点，论坛里已经有不少相关的帖子和程序，我想做的是反距离权重法插值，格点资料是2.5x2.5的。大家可能都有这样的体会：别人的程序看不懂……之前下载的双线性插值法程序也是因为不明白（记不住）所有变量的含义所以看不下去。理都懂，所以想试试自己编一个看看。
我自己因为本科学过FORTRAN，并且大二、三、四、研，学校设置的课程、实习以及老师的安排都需要用FORTRAN，所以没有学完就丢掉，还是比较庆幸的。基础不算太差，但是我觉得大部分问题都能只用do循环、if语句等解决。

下面讲我编程的拙见：
首先确定问题。下面是在草稿本上写的，画了坐标示意图，公式很明确，要确定四个格点是主要的问题，尝试将四个格点的经纬度分别减去站点经纬度，取绝对值后，小于2.5，认为是离站点最近的四个格点。思路要清晰，记得以前老师上课时说写程序之前要画流程图，虽然没有做到那么仔细，但是一上手就写程序肯定是很困难的。

然后尝试写主要的步骤，格点乘上2.5，就是经纬度的值，与站点经纬度相减并取绝对值，简单的do和if语句就能搞定

之后就试着运行，试着输入台站号和经纬度，运行出了点问题，但是四个点都找到了，这就是最主要的部分，剩下的部分就是按公式，四则运算，修改修改修改，更好写了。

这是花了半个小时左右的成果，比起看不懂程序一天天不想做，来得更有动力些，希望能鼓励鼓励大家！
个人能力有限，一些见解不到位或者程序有错误请大力指正！！！
此程序！！！***局限于单个站点多时次插值***！！！！
2017年1月3号：收到论坛朋友的指正，发现了一些公式以及计算方法的问题，特此更正！

implicit none
real re(144,73,67),v1(67)
!re(144,73,67)是经过计算的再分析资料67年年平均值（1984-2014），2.5*2.5,纬度-90~90，经度0~360
!后面有说明负值纬度的修正，v1(67)是即将计算的新插值序列
real d,lat,lon,fenzi(67),fenmu
!d是距离的平方，因为计算距离d=sqrt(((i-1)*2.5-lon)**2+((j-1)*2.5-lat)**2),之后公式还是有平方，干脆不开方也不平方。
!lat是输入的站点纬度，lon是输入的站点精度。此程序！！！***局限于单个站点多时次插值***！！！！
!fenzi(67)是插值公式的分子，fenmu是插值公式的分母
integer i,j,t
open(10,file='d:\wind\201602\data\1948_2014yrmean.grd',form='binary')
open(20,file='d:\wind\201602\data\beijingchazhi.txt') !北京插值为例
read(*,*) lon,lat
do t=1,67
read(10)((re(i,j,t),i=1,144),j=1,73)
enddo
fenmu=0
do j=1,73
do i=1,144
if(abs((i-1)*2.5-lon)<=2.5.and.abs((j-1)*2.5-lat)<=2.5)then
!判断与站点经纬度差值的绝对值小于2.5视为最近的四个点。这里-1的原因是，计算格点
!的公式是，以2.5度、0~360经度为例，(360-0)/2.5+1=144。所以-1体现在此。
d=1/((((i-1)*2.5-lon)*111)**2+(((j-1)*2.5-lat)*111)**2) !距离的平方的倒数，同时格点1度大约是111千米
print*,d
fenmu=fenmu+d
do t=1,67
fenzi(t)=fenzi(t)+d*re(i,j+36,t)
!j+36是将负值纬度调到正值。北京39.7N为例，计算得到相近两个纬向格点是16、17，但这前提是以赤道为格点1的起点，
!而-90~0度格点为1-37，因此要以格点37为赤道起点的话，16、17就要加上1和37的差值36。我说得不清楚或者不明白可以
!用grads自己试试，set x 16得到的纬度是-52.5，而加上36后，set x 52得到纬度37.5
enddo
endif
enddo
enddo
do t=1,67
v1(t)=fenzi(t)/fenmu
enddo
write(20,'(1x,f10.5)') (v1(t),t=1,67)
end

复制代码

ishitheno1 · 发表于 2016-1-16 19:52:27

此处必须有掌声！加油

subtropical · 发表于 2016-1-16 20:14:11

{:5_213:}{:5_213:}为坚持不懈的人们鼓掌

四叶草 · 发表于 2016-1-16 20:17:49

必须鼓掌

hwl1943 · 发表于 2016-1-16 20:50:41

给楼主掌声，加油

topmad · 发表于 2016-1-16 20:55:25

赞一个
楼主的方式我很喜欢先数学化然后在纸上程序化最后再电脑里一下基本上就好了
写上注释一切都好说

andrewsoong · 发表于 2016-1-16 21:12:42

超赞！！！鼓励！！！

自导自演 · 发表于 2016-1-16 22:14:43

楼主加油，棒

苏叶绿阳 · 发表于 2016-1-17 17:23:08

鼓励鼓励，

咕噜 · 发表于 2016-1-17 21:03:31

申请贡献呀~@mofangbao

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

[经验总结] 分享我写FORTRAN程序（反距离权重法插值）的一些步骤和心得

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

评分

本帖被以下淘专辑推荐: