请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 2115|回复: 5

[混合编程] IDL之通用版克里金Kriging插值程序

[复制链接]

新浪微博达人勋

发表于 2020-7-2 17:23:18 | 显示全部楼层 |阅读模式

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

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 15195775117 于 2020-12-24 09:39 编辑

飞鸽截图20190619141710.png


一、缘起

克里金法是我工作之后才知道的插值法,其插值效果极其优美,是我平时用的最多的插值法。

二、特性

克里金法(Kriging)是依据协方差函数随机过程/随机场进行空间建模和预测(插值)的回归算法。在特定的随机过程,例如固有平稳过程中,克里金法能够给出最优线性无偏估计,因此在地统计学中也被称为空间最优无偏估计器,被应用于地理科学环境科学大气科学等领域。
克里金法是典型的地统计学算法,以鄙人浅显的经验,它的插值结果与原始散点的值并不是绝对对应,而是整体一致但可能有局部差别

三、IDL中的克里金

IDL中有二维克里金插值函数Krig2D,但并不是全自动,克里金插值法有几种常用插值模型,需要人为提供模型参数。

我也买了2本专业书来研究克里金,但是推导公式太麻烦了,于是我习惯于使用help中提供的参数。
没有了配置模型参数的烦恼,克里金法的参数只有经度、纬度、值3者,根据经验,一个位置不能有多个值,于是我写了去除重复位置的函数,对于同一个位置的值,做取平均处理。

四、去除重复位置的函数

;通用Krig的辅助函数
;剔除重复[x,y]位置的程序,lat,lon,ext是向量
;对于同一的位置的值,取平均
function lat_lon_ext,lat,lon,ext
  ;如果向量大小不一致就报错:
  n1=n_elements(ext)
  n2=n_elements(lon)
  n3=n_elements(lat)
  if(n1 ne n2 or n2 ne n3)then begin
    warning=dialog_message('数量不一致!')
    stop
  endif

  n=n1;数量一致用n表示
  book=[-123,-123,-123,-123]
  for i=0,n-1 do begin
    p=where(lon eq lon and lat eq lat,count)
    book=[[book],[lat,lon,count,mean(ext[p])]]
  endfor
  book=book[*,1:-1]
  sizebook=size(book)
  book2=[-123,-123,-123,-123]
  for i=0,sizebook[2]-1 do begin
    p=where(book2[0,*] eq book[0,i] and book2[1,*] eq book[1,i],count)
    if(count eq 0)then begin
      book2=[[book2],[book[*,i]]]
    endif
    if count eq 1 then continue
    if(count gt 1 or count lt 0)then begin
      warning=DIALOG_MESSAGE('出错!')
      stop
    endif
  endfor
  book2=book2[*,1:-1]
  lon2=transpose(book2[1,*])
  lat2=transpose(book2[0,*])
  value2=transpose(book2[3,*])
  re={lon:lon2,lat:lat2,value:value2}
  return,re
end

五、通用克里金函数

;返回值结构体中增设了fail用于判断插值是否失败
function GeneralKrig,lon1,lat1,value1
  ;插值成败指示变量fail:
  fail=0
  ;去除重复地点值:
  k=lat_lon_ext(lat1,lon1,value1)
  ;krig插值,使用了IDL的krig2d的help中的模型参数:
  value2=KRIG2D(k.value,k.lon,k.lat, EXPONENTIAL=[0.5,0.2,1], $
    NX=200, NY=300, XOUT=xout, YOUT=yout);NX和NY表示返回值数组的大小
  ;如果value2中有nan值,则插值失败:
  p=where(finite(value2) eq 0,count)
  if(count gt 0)then begin
    warning=dialog_message('克里金插值失败!')
    fail=1
  endif
  re={lon:xout,lat:yout,value:value2,fail:fail}
  return,re
end

六、调用方法

去除重复点的lat_lon_ext函数不用管,它是给GeneralKrig使用的,
re=GeneralKrig(经度序列,纬度序列,值序列)
对于返回值先看re.fail是否是0,0是成功,1是插值失败
插值后的经度数组是re.lon
插值后的纬度数组是re.lat
插值后的值数组是re.value
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2021-10-9 16:53:04 | 显示全部楼层
您真是好厉害!这也是我要学的!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2021-10-10 14:59:43 | 显示全部楼层
赵家清欢 发表于 2021-10-9 16:53
您真是好厉害!这也是我要学的!

过誉了,共同进步!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2022-2-12 13:10:35 | 显示全部楼层
共同进步!
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2022-7-17 08:36:39 | 显示全部楼层
EXPONENTIAL如何调整,有没有推荐呢?
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2022-7-17 10:58:05 | 显示全部楼层
金今 发表于 2022-7-17 08:36
EXPONENTIAL如何调整,有没有推荐呢?

这里的模型参数我也不太懂
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

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

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

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