爱气象，爱气象家园！ 
·设为首页 开启辅助访问 充值积分新帖子新回复热门贴精华帖

气象家园

 找回密码
 立即注册

QQ登录

只需一步，快速开始

搜索
气象家园»论坛 › 气象软件开发 › c#/asp.net › C#编写格点插值法
1234下一页
返回列表 发新帖
查看: 16899|回复: 36

C#编写格点插值法

[复制链接]
tmmiw
发表于 2014-8-5 08:57:08 | 显示全部楼层 |阅读模式

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

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

x
为了步浪费您的金钱直接附上源码.注意本源码用C#编写.
   static int NS = 554, MX = 30, MY = 24, ND = 1;
        //NS离散点数目，观测站总数，MX精度格点数，MY纬度格点数
        float[] lt = new float[NS];
        float[] ln = new float[NS];
        float[] lat = new float[MY];
        float[] lon = new float[MX];
        float[,] p1x1 = new float[MX,MY];
        float[,] t1x1 = new float[MX, MY];
        float[,] ipr = new float[NS, ND];
        float[,] t = new float[NS, ND];

            int nlev, nlfag;
        float tim=Convert.ToSingle(0.25);//插值半径0.25度




       lon[0] = Convert.ToSingle(99.8);//经度开始值99.8-107.5
            lat[0] = Convert.ToSingle(26.5);//纬度开始值26.5-32.65
            for (int i = 1; i < MX; i++)
                lon[i] = lon[i - 1] + Convert.ToSingle(1.0);//给格点上的经度负值
            for (int i = 1; i < MY; i++)
                lat[i] = lat[i - 1] + Convert.ToSingle(1.0);//给格点上的纬度赋值

  //初始化p1x1
            for (int j = 0; j < MY; j++)
                for (int i = 0; i < MX; i++)
                    p1x1[i, j] = Convert.ToSingle(-999.9);
            //开始计算了

            for (int it = 0; it < ND; it++)
            {
                #region 格点插值
                float wt = Convert.ToSingle(0.0), r1x1 = Convert.ToSingle(0.0), w = Convert.ToSingle(0.0);
                for (int j = 1; j < MY - 1; j++)
                {
                    for (int i = 1; i < MX - 1; i++)
                    {
                        float x0 = lon[i];
                        float y0 = lat[j];
                        wt = Convert.ToSingle(0.0); r1x1 = Convert.ToSingle(0.0); w = Convert.ToSingle(0.0);
                        for (int ii = 0; ii < nlev; ii++)
                        {
                            if (ipr[ii, it] != -999.9)
                            {
                                float x = ln[ii];
                                float y = lt[ii];
                                float d = Convert.ToSingle(Math.Sqrt((x - x0) * (x - x0) + (y - y0) * (y - y0)));
                                if (d <1.25* tim)
                                {
                                    if (d < 0.1*tim)
                                        w = Convert.ToSingle(100.0);//避免间隔太近的点贡献过大
                                    else
                                        w = Convert.ToSingle(1.0) / (d * d);
                                    r1x1 = r1x1 + ipr[ii, it] * w;
                                    wt += w;
                                }
                            }
                        }
                        if (wt < 0.2)
                            p1x1[i, j] = Convert.ToSingle(-999.9);
                        else
                            p1x1[i, j] = r1x1 / wt;

                    }
                }
            #endregion

//离散格点插值就完成了,但在一些站点稀疏的地方可能有的格点还没有附上值
//下面开始用格点插值格点
  #region 跨一个格点补充缺数点
                for (int j = 1; j < MY - 1; j++)
                {
                    for (int i = 1; i < MX - 1; i++)
                    {
                        if(p1x1[i,j]==-999.99 &&
                            (p1x1[i+1,j]!=-999.9 || p1x1[i+1,j+1]!=-999.9 || p1x1[i+1,j-1]!=-999.9) &&
                            (p1x1[i-1,j]!=-999.9 || p1x1[i-1,j+1]!=-999.9 || p1x1[i-1,j-1]!=-999.9) &&
                            (p1x1[i,j+1]!=-999.9 || p1x1[i+1,j+1]!=-999.9 || p1x1[i-1,j+1]!=-999.9) &&
                            (p1x1[i,j-1]!=-999.9 || p1x1[i+1,j-1]!=-999.9 || p1x1[i-1,j-1]!=-999.9) &&
                            (p1x1[i+1,j+1]!=-999.9 || p1x1[i+1,j]!=-999.9 || p1x1[i,j+1]!=-999.9) &&
                            (p1x1[i-1,j-1]!=-999.9 || p1x1[i-1,j]!=-999.9 || p1x1[i,j-1]!=-999.9) &&
                            (p1x1[i+1,j-1]!=-999.9 || p1x1[i+1,j]!=-999.9 || p1x1[i,j-1]!=-999.9) &&
                            (p1x1[i-1,j+1]!=-999.9 || p1x1[i-1,j]!=-999.9 || p1x1[i,j+1]!=-999.9))
                        {
                              wt = Convert.ToSingle(0.0);r1x1 = Convert.ToSingle(0.0); w = Convert.ToSingle(0.0);
                            for(int n=j-1;n<j+2;n++)
                            {
                                for(int m=i-1;m<i+2;n++)
                                {
                                    if(p1x1[m,n]!=-999.9)
                                    {
                                        w=Convert.ToSingle(1.0/Convert.ToSingle((i-m)*(i-m)+(j-n)*(j-n)));
                                        wt=wt+w;
                                        r1x1=r1x1 +p1x1[m,n]*w;

                                    }
                                }
                            }
                            p1x1[i,j]=r1x1/wt;
                        }
                    }
                }
                #endregion

   #region 跨二个格点补充缺数点
                for (int j = 2; j < MY - 2; j++)
                {
                    for (int i = 2; i < MX - 2; i++)
                    {
                        if(p1x1[i,j]==-999.99 &&
                            (p1x1[i+2,j]!=-999.9 || p1x1[i+2,j+2]!=-999.9 || p1x1[i+2,j-2]!=-999.9) &&
                            (p1x1[i-2,j]!=-999.9 || p1x1[i-2,j+2]!=-999.9 || p1x1[i-2,j-2]!=-999.9) &&
                            (p1x1[i,j+2]!=-999.9 || p1x1[i+2,j+2]!=-999.9 || p1x1[i-2,j+2]!=-999.9) &&
                            (p1x1[i,j-2]!=-999.9 || p1x1[i+2,j-2]!=-999.9 || p1x1[i-2,j-2]!=-999.9) &&
                            (p1x1[i+2,j+2]!=-999.9 || p1x1[i+2,j]!=-999.9 || p1x1[i,j+2]!=-999.9) &&
                            (p1x1[i-2,j-2]!=-999.9 || p1x1[i-2,j]!=-999.9 || p1x1[i,j-2]!=-999.9) &&
                            (p1x1[i+2,j-2]!=-999.9 || p1x1[i+2,j]!=-999.9 || p1x1[i,j-2]!=-999.9) &&
                            (p1x1[i-2,j+2]!=-999.9 || p1x1[i-2,j]!=-999.9 || p1x1[i,j+2]!=-999.9))
                        {
                              wt = Convert.ToSingle(0.0); r1x1 = Convert.ToSingle(0.0); w = Convert.ToSingle(0.0);
                            for(int n=j-2;n<j+3;n++)
                            {
                                for(int m=i-2;m<i+3;n++)
                                {
                                    if(p1x1[m,n]!=-999.9)
                                    {
                                        w = Convert.ToSingle(1.0 / Convert.ToSingle((i - m) * (i - m) + (j - n) * (j - n)));
                                        wt=wt+w;
                                        r1x1=r1x1 +p1x1[m,n]*w;

                                    }
                                }
                            }
                            p1x1[i,j]=r1x1/wt;
                        }
                    }
                }
                #endregion

#region 跨三个格点补充缺数点
                for (int j =3; j < MY -3; j++)
              {
                  for (int i =3; i < MX -3; i++)
                  {
                      if (p1x1[i, j] == -999.99 &&
                          (p1x1[i +3, j] != -999.9 || p1x1[i +3, j +3] != -999.9 || p1x1[i +3, j -3] != -999.9) &&
                          (p1x1[i -3, j] != -999.9 || p1x1[i -3, j +3] != -999.9 || p1x1[i -3, j -3] != -999.9) &&
                          (p1x1[i, j +3] != -999.9 || p1x1[i +3, j +3] != -999.9 || p1x1[i -3, j +3] != -999.9) &&
                          (p1x1[i, j -3] != -999.9 || p1x1[i +3, j -3] != -999.9 || p1x1[i -3, j -3] != -999.9) &&
                          (p1x1[i +3, j +3] != -999.9 || p1x1[i +3, j] != -999.9 || p1x1[i, j +3] != -999.9) &&
                          (p1x1[i -3, j -3] != -999.9 || p1x1[i -3, j] != -999.9 || p1x1[i, j -3] != -999.9) &&
                          (p1x1[i +3, j -3] != -999.9 || p1x1[i +3, j] != -999.9 || p1x1[i, j -3] != -999.9) &&
                          (p1x1[i -3, j +3] != -999.9 || p1x1[i -3, j] != -999.9 || p1x1[i, j +3] != -999.9))
                      {
                          wt = Convert.ToSingle(0.0); r1x1 = Convert.ToSingle(0.0); w = Convert.ToSingle(0.0);
                          for (int n = j -3; n < j + 4; n++)
                          {
                              for (int m = i -3; m < i + 4; n++)
                              {
                                  if (p1x1[m, n] != -999.9)
                                  {
                                      w = Convert.ToSingle(1.0 / Convert.ToSingle((i - m) * (i - m) + (j - n) * (j - n)));
                                      wt = wt + w;
                                      r1x1 = r1x1 + p1x1[m, n] * w;

                                  }
                              }
                          }
                          p1x1[i, j] = r1x1 / wt;
                      }
                  }
              }
                #endregion

//再做一遍跨两个和一个格点补充
#region 跨二个格点补充缺数点
                for (int j = 2; j < MY - 2; j++)
                {
                    for (int i = 2; i < MX - 2; i++)
                    {
                        if (p1x1[i, j] == -999.99 &&
                            (p1x1[i + 2, j] != -999.9 || p1x1[i + 2, j + 2] != -999.9 || p1x1[i + 2, j - 2] != -999.9) &&
                            (p1x1[i - 2, j] != -999.9 || p1x1[i - 2, j + 2] != -999.9 || p1x1[i - 2, j - 2] != -999.9) &&
                            (p1x1[i, j + 2] != -999.9 || p1x1[i + 2, j + 2] != -999.9 || p1x1[i - 2, j + 2] != -999.9) &&
                            (p1x1[i, j - 2] != -999.9 || p1x1[i + 2, j - 2] != -999.9 || p1x1[i - 2, j - 2] != -999.9) &&
                            (p1x1[i + 2, j + 2] != -999.9 || p1x1[i + 2, j] != -999.9 || p1x1[i, j + 2] != -999.9) &&
                            (p1x1[i - 2, j - 2] != -999.9 || p1x1[i - 2, j] != -999.9 || p1x1[i, j - 2] != -999.9) &&
                            (p1x1[i + 2, j - 2] != -999.9 || p1x1[i + 2, j] != -999.9 || p1x1[i, j - 2] != -999.9) &&
                            (p1x1[i - 2, j + 2] != -999.9 || p1x1[i - 2, j] != -999.9 || p1x1[i, j + 2] != -999.9))
                        {
                            wt = Convert.ToSingle(0.0); r1x1 = Convert.ToSingle(0.0); w = Convert.ToSingle(0.0);
                            for (int n = j - 2; n < j + 3; n++)
                            {
                                for (int m = i - 2; m < i + 3; n++)
                                {
                                    if (p1x1[m, n] != -999.9)
                                    {
                                        w = Convert.ToSingle(1.0 / Convert.ToSingle((i - m) * (i - m) + (j - n) * (j - n)));
                                        wt = wt + w;
                                        r1x1 = r1x1 + p1x1[m, n] * w;

                                    }
                                }
                            }
                            p1x1[i, j] = r1x1 / wt;
                        }
                    }
                }
                #endregion

#region 跨一个格点补充缺数点
                for (int j = 1; j < MY - 1; j++)
                {
                    for (int i = 1; i < MX - 1; i++)
                    {
                        if (p1x1[i, j] == -999.99 &&
                            (p1x1[i + 1, j] != -999.9 || p1x1[i + 1, j + 1] != -999.9 || p1x1[i + 1, j - 1] != -999.9) &&
                            (p1x1[i - 1, j] != -999.9 || p1x1[i - 1, j + 1] != -999.9 || p1x1[i - 1, j - 1] != -999.9) &&
                            (p1x1[i, j + 1] != -999.9 || p1x1[i + 1, j + 1] != -999.9 || p1x1[i - 1, j + 1] != -999.9) &&
                            (p1x1[i, j - 1] != -999.9 || p1x1[i + 1, j - 1] != -999.9 || p1x1[i - 1, j - 1] != -999.9) &&
                            (p1x1[i + 1, j + 1] != -999.9 || p1x1[i + 1, j] != -999.9 || p1x1[i, j + 1] != -999.9) &&
                            (p1x1[i - 1, j - 1] != -999.9 || p1x1[i - 1, j] != -999.9 || p1x1[i, j - 1] != -999.9) &&
                            (p1x1[i + 1, j - 1] != -999.9 || p1x1[i + 1, j] != -999.9 || p1x1[i, j - 1] != -999.9) &&
                            (p1x1[i - 1, j + 1] != -999.9 || p1x1[i - 1, j] != -999.9 || p1x1[i, j + 1] != -999.9))
                        {
                            wt = Convert.ToSingle(0.0); r1x1 = Convert.ToSingle(0.0); w = Convert.ToSingle(0.0);
                            for (int n = j - 1; n < j + 2; n++)
                            {
                                for (int m = i - 1; m < i + 2; n++)
                                {
                                    if (p1x1[m, n] != -999.9)
                                    {
                                        w = Convert.ToSingle(1.0 / Convert.ToSingle((i - m) * (i - m) + (j - n) * (j - n)));
                                        wt = wt + w;
                                        r1x1 = r1x1 + p1x1[m, n] * w;

                                    }
                                }
                            }
                            p1x1[i, j] = r1x1 / wt;
                        }
                    }
                }
                #endregion

//程序到这里就结束了,这样一方面可以避免格点缺值,也消除过于平滑看不到峰值和谷值得现象.

评分

参与人数 1金钱 +12 贡献 +3 收起 理由
mofangbao + 12 + 3

查看全部评分

密码修改失败请联系微信：mofangbao
回复

使用道具 举报

mofangbao
0
早起挑战累计收入
发表于 2014-8-5 13:18:48 | 显示全部楼层
楼主介绍下用的是啥算法？还是自己原创的算法呢
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

llsshh1985
发表于 2014-8-5 14:07:51 | 显示全部楼层
如果是原创的，真是厉害啊
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

tmmiw
 楼主| 发表于 2014-8-5 14:10:55 | 显示全部楼层
程序是自己写的,算法当然不是,就是离散点格点化.
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

tmmiw
 楼主| 发表于 2014-8-5 14:39:14 | 显示全部楼层
原理请参看：

格点数据插值(算法)
http://bbs.06climate.com/forum.p ... 3&fromuid=37658
(出处: 气象家园)
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

fengyu1983
发表于 2014-8-20 17:38:30 | 显示全部楼层
这论坛上强者太多了
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

y18735399
发表于 2014-12-26 18:07:50 | 显示全部楼层
看别人的代码是很痛苦的事情
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

liao001
发表于 2015-1-3 17:40:20 | 显示全部楼层
程序跑起来么？
nlev 都没有给值。W算不出来哦。
ND是什么？
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

liao001
发表于 2015-1-3 17:50:30 | 显示全部楼层
能贴出完整能运行的源码来看一下么？
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

虎皮
发表于 2015-1-18 09:29:06 | 显示全部楼层
楼主的代码可读性太差啊，错误条件先返回在哪里？-999.9写了无数次你不烦？
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

下一页 »
1234下一页
返回列表 发新帖
高级模式
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

小黑屋|帮助中心|清除痕迹

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表，不代表气象家园立场，禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表