python 普通克里金（Kriging）法的实现

克里金法时一种用于空间插值的地学统计方法。

克里金法用半变异测定空间要素，要素即自相关要素。

半变异公式为：

其中γ(h) 是已知点 xi 和 xj 的半变异，***h***表示这两个点之间的距离，z是属性值。

假设不存在漂移，普通克里金法重点考虑空间相关因素，并用拟合的半变异直接进行插值。

估算某测量点z值的通用方程为：

式中，z0是待估计值，zx是已知点x的值，Wx是每个已知点关联的权重，s是用于估计的已知点数目。
权重可以由一组矩阵方程得到。

此程序对半变异进行拟合时采用的时最简单的正比例函数拟合

数据为csv格式

保存格式如下：

第一行为第一个点以此类推

最后一行是待求点坐标，其中z为未知值，暂且假设为0

代码如下：

import numpy as np
from math import*
from numpy.linalg import *
h_data=np.loadtxt(open('高程点数据.csv'),delimiter=",",skiprows=0)
print('原始数据如下(x,y,z):\n未知点高程初值设为0\n',h_data)
def dis(p1,p2):
 a=pow((pow((p1[0]-p2[0]),2)+pow((p1[1]-p2[1]),2)),0.5)
 return a
def rh(z1,z2):
 r=1/2*pow((z1[2]-z2[2]),2)
 return r
def proportional(x,y):
 xx,xy=0,0
 for i in range(len(x)):
  xx+=pow(x[i],2)
  xy+=x[i]*y[i]
 k=xy/xx
 return k
r=[];pp=[];p=[];
for i in range(len(h_data)):
 pp.append(h_data[i])
for i in range(len(pp)):
 for j in range(len(pp)):
  p.append(dis(pp[i],pp[j]))
  r.append(rh(pp[i],pp[j]))
r=np.array(r).reshape(len(h_data),len(h_data))
r=np.delete(r,len(h_data)-1,axis =0)
r=np.delete(r,len(h_data)-1,axis =1)

h=np.array(p).reshape(len(h_data),len(h_data))
h=np.delete(h,len(h_data)-1,axis =0)
oh=h[:,len(h_data)-1]
h=np.delete(h,len(h_data)-1,axis =1)

hh=np.triu(h,0)
rr=np.triu(r,0)
r0=[];h0=[];
for i in range(len(h_data)-1):
 for j in range(len(h_data)-1):
  if hh[i][j] !=0:
   a=h[i][j]
   h0.append(a)
  if rr[i][j] !=0:
   a=rr[i][j]
   r0.append(a)
k=proportional(h0,r0)
hnew=h*k
a2=np.ones((1,len(h_data)-1))
a1=np.ones((len(h_data)-1,1))
a1=np.r_[a1,[[0]]]
hnew=np.r_[hnew,a2]
hnew=np.c_[hnew,a1]
print('半方差联立矩阵：\n',hnew)
oh=np.array(k*oh)
oh=np.r_[oh,[1]]
w=np.dot(inv(hnew),oh)
print('权阵运算结果：\n',w)
z0,s2=0,0
for i in range(len(h_data)-1):
 z0=w[i]*h_data[i][2]+z0
 s2=w[i]*oh[i]+s2
s2=s2+w[len(h_data)-1]
print('未知点高程值为：\n',z0)
print('半变异值为：\n',pow(s2,0.5))
input()

运算结果

python初学，为了完成作业写了个小程序来帮助计算，因为初学知识有限，有很多地方写的很复杂，可以优化的地方很多。 还望读者谅解，欢迎斧正谢谢！

参考文献：
【1】(美)张康聪 著;陈健飞等译. 地理信息系统导论（第三版）. 北京：清华大学出版社, 2009.04.

以上就是本文的全部内容，希望对大家的学习有所帮助，也希望大家多多支持我们。

赞 (0)

相关推荐

  • python编写的最短路径算法

    python编写的最短路径算法

    一心想学习算法,很少去真正静下心来去研究,前几天趁着周末去了解了最短路径的资料,用python写了一个最短路径算法.算法是基于带权无向图去寻找两个点之间的最短路径,数据存储用邻接矩阵记录.首先画出一幅无向图如下,标出各个节点之间的权值. 其中对应索引: A --> 0 B--> 1 C--> 2 D-->3 E--> 4 F--> 5 G--> 6 邻接矩阵表示无向图: 算法思想是通过Dijkstra算法结合自身想法实现的.大致思路是:从起始点开始,搜索周围的路径

  • python使用rsa加密算法模块模拟新浪微博登录

    PC登录新浪微博时,在客户端用js预先对用户名.密码都进行了加密,而且在POST之前会GET一组参数,这也将作为POST_DATA的一部分.这样,就不能用通常的那种简单方法来模拟POST登录(比如人人网). 通过爬虫获取新浪微博数据,模拟登录是必不可少的. 1.在提交POST请求之前,需要GET获取四个参数(servertime,nonce,pubkey和rsakv),不是之前提到的只是获取简单的servertime,nonce,这里主要是由于js对用户名.密码加密方式改变了. 1.1 由于加密

  • Python实现的凯撒密码算法示例

    Python实现的凯撒密码算法示例

    本文实例讲述了Python实现的凯撒密码算法.分享给大家供大家参考,具体如下: 一 介绍 凯撒密码是一种非常古老的加密方法,相传当年凯撒大地行军打仗时为了保证自己的命令不被敌军知道,就使用这种特殊的方法进行通信,以确保信息传递的安全.他的原理很简单,说到底就是字母于字母之间的替换.下面让我们看一个简单的例子:"baidu"用凯撒密码法加密后字符串变为"edlgx",它的原理是什么呢?把"baidu"中的每一个字母按字母表顺序向后移3位,所得的结果

  • python k-近邻算法实例分享

    简单说明 这个算法主要工作是测量不同特征值之间的距离,有个这个距离,就可以进行分类了. 简称kNN. 已知:训练集,以及每个训练集的标签. 接下来:和训练集中的数据对比,计算最相似的k个距离.选择相似数据中最多的那个分类.作为新数据的分类. python实例 复制代码 代码如下: # -*- coding: cp936 -*- #win系统中应用cp936编码,linux中最好还是utf-8比较好.from numpy import *#引入科学计算包import operator #经典pyt

  • 手把手教你python实现SVM算法

    手把手教你python实现SVM算法

    什么是机器学习 (Machine Learning) 机器学习是研究计算机怎样模拟或实现人类的学习行为,以获取新的知识或技能,重新组织已有的知识结构使之不断改善自身的性能.它是人工智能的核心,是使计算机具有智能的根本途径,其应用遍及人工智能的各个领域. 机器学习的大致分类: 1)分类(模式识别):要求系统依据已知的分类知识对输入的未知模式(该模式的描述)作分析,以确定输入模式的类属,例如手写识别(识别是不是这个数). 2)问题求解:要求对于给定的目标状态,寻找一个将当前状态转换为目标状态的动作序

  • python实现协同过滤推荐算法完整代码示例

    python实现协同过滤推荐算法完整代码示例

    测试数据 http://grouplens.org/datasets/movielens/ 协同过滤推荐算法主要分为: 1.基于用户.根据相邻用户,预测当前用户没有偏好的未涉及物品,计算得到一个排序的物品列表进行推荐 2.基于物品.如喜欢物品A的用户都喜欢物品C,那么可以知道物品A与物品C的相似度很高,而用户C喜欢物品A,那么可以推断出用户C也可能喜欢物品C. 不同的数据.不同的程序猿写出的协同过滤推荐算法不同,但其核心是一致的: 1.收集用户的偏好 1)不同行为分组 2)不同分组进行加权计算用

  • python实现RSA加密(解密)算法

    python实现RSA加密(解密)算法

    RSA是目前最有影响力的公钥加密算法,它能够抵抗到目前为止已知的绝大多数密码攻击,已被ISO推荐为公钥数据加密标准. 今天只有短的RSA钥匙才可能被强力方式解破.到2008年为止,世界上还没有任何可靠的攻击RSA算法的方式.只要其密钥的长度足够长,用RSA加密的信息实际上是不能被解破的.但在分布式计算和量子计算机理论日趋成熟的今天,RSA加密安全性受到了挑战. RSA算法基于一个十分简单的数论事实:将两个大素数相乘十分容易,但是想要对其乘积进行因式分解却极其困难,因此可以将乘积公开作为加密密钥.

  • python冒泡排序算法的实现代码

    1.算法描述:(1)共循环 n-1 次(2)每次循环中,如果 前面的数大于后面的数,就交换(3)设置一个标签,如果上次没有交换,就说明这个是已经好了的. 2.python冒泡排序代码 复制代码 代码如下: #!/usr/bin/python# -*- coding: utf-8 -*- def bubble(l):    flag = True    for i in range(len(l)-1, 0, -1):        if flag:             flag = False

  • python实现识别手写数字 python图像识别算法

    python实现识别手写数字 python图像识别算法

    写在前面 这一段的内容可以说是最难的一部分之一了,因为是识别图像,所以涉及到的算法会相比之前的来说比较困难,所以我尽量会讲得清楚一点. 而且因为在编写的过程中,把前面的一些逻辑也修改了一些,将其变得更完善了,所以一切以本篇的为准.当然,如果想要直接看代码,代码全部放在我的GitHub中,所以这篇文章主要负责讲解,如需代码请自行前往GitHub. 本次大纲 上一次写到了数据库的建立,我们能够实时的将更新的训练图片存入CSV文件中.所以这次继续往下走,该轮到识别图片的内容了. 首先我们需要从文件夹中

  • python 普通克里金（Kriging）法的实现

    克里金法时一种用于空间插值的地学统计方法. 克里金法用半变异测定空间要素,要素即自相关要素. 半变异公式为: 其中γ(h) 是已知点 xi 和 xj 的半变异,***h***表示这两个点之间的距离,z是属性值. 假设不存在漂移,普通克里金法重点考虑空间相关因素,并用拟合的半变异直接进行插值. 估算某测量点z值的通用方程为: 式中,z0是待估计值,zx是已知点x的值,Wx是每个已知点关联的权重,s是用于估计的已知点数目. 权重可以由一组矩阵方程得到. 此程序对半变异进行拟合时采用的时最简单的正比例

  • python清除字符串里非字母字符的方法

    本文实例讲述了python清除字符串里非字母字符的方法.分享给大家供大家参考.具体如下: s = "hello world! how are you? 0" # Short version print filter(lambda c: c.isalpha(), s) # Faster version for long ASCII strings: id_tab = "".join(map(chr, xrange(256))) tostrip = "&quo

  • python清除字符串里非数字字符的方法

    本文实例讲述了python清除字符串里非数字字符的方法.分享给大家供大家参考.具体如下: import re s = "how19 a*re 254y**ou?" # Using regular expressions print re.sub("\D", "", s) 希望本文所述对大家的Python程序设计有所帮助.

  • python将MongoDB里的ObjectId转换为时间戳的方法

    本文实例讲述了python将MongoDB里的ObjectId转换为时间戳的方法.分享给大家供大家参考.具体分析如下: MongoDB里的_id字段前四位是时间戳的16进制表示,通过Python可以很容易从_id中提取出时间戳来 def timestamp_from_objectid(objectid): result = 0 try: result = time.mktime(objectid.generation_time.timetuple()) except: pass return r

  • Python 在OpenCV里实现仿射变换—坐标变换效果

    在现实的图像操作软件中,经常碰到的不是给出放大多少倍,而是由用户在软件的界面上选择多大的区域,或者选择几个点,那么这样情况下,怎么样来计算出变换矩阵呢?从前面知道变换矩阵是2X3的矩阵,说明有六个未知数,又有中学的代数知识知道要解决六个未知数,那么方程组至少要联立三条方程,要准备三条方程的先决条件,就是要有三组坐标.因此,只要在用户选择的区域里找到三个不同点的坐标,就可以计算出变换矩阵.如果给出三组坐标[0, 0], [200, 0], [0, 200],通过变换之后新坐标是[0, 0], [1

  • python在OpenCV里实现投影变换效果

    前面学习了仿射变换,是经常使用到的变换,也很容易理解.在日常生活中,经常会遇到下面这种的情况: 仔细地观察比亚迪秦这台汽车的车牌,发现它拍照的角度不是垂直的方向,而是有一个角度,当要进行车牌识别的时候,发现字符是变形的,与电脑里比较的图片肯定有区别,因此识别不出来.这时怎么办呢?就需要经过一个投影变换才可以把车牌号纠正过来,才能进入识别过程. 好吧,到这里认识到投影变换的感性认识了,那么你又会继续考虑下一个问题,在软件里怎么样计算呢,难道还是使用仿射变换的矩阵.从这里看一下,前面闽A比较大,后面

  • Python在OpenCV里实现极坐标变换功能

    在中学里学习过直角坐标系,也叫做笛卡尔坐标系,它是正交坐标系,不过也学习过极坐标系,这种坐标系比较适合大炮发射的场合.极坐标系的定义如下: 在 平面内取一个定点O, 叫极点,引一条射线Ox,叫做极轴,再选定一个长度单位和角度的正方向(通常取逆时针方向).对于平面内任何一点M,用ρ表示线段OM的长度,θ表示从Ox到OM的角度,ρ叫做点M的极径,θ叫做点M的极角,有序数对 (ρ,θ)就叫点M的极坐标,这样建立的坐标系叫做极坐标系. 极坐标很方便应用到雷达上面,因为雷达不断地转动,反射回来的波计算出距

  • 解决python给列表里添加字典时被最后一个覆盖的问题

    如下所示: >>> item={} ; items=[] #先声明一个字典和一个列表,字典用来添加到列表里面 >>> item['index']=1 #给字典赋值 >>> items.append(item) >>> items [{'index': 1}] #添加到列表里面复合预期 >>> item['index']=2 #现在修改字典 >>> item {'index': 2} #修改成功 &g

  • 解决Python中list里的中文输出到html模板里的问题

    最仅在做一个数据分析的功能时候遇到将list中的中文字符按照数组的形式输出到html模板里的js中进行处理,但是直接输出模板会按照unicode编码输出,这个问题真的让人头大. 本方法实在flask框架里完成的,以下是解决方法,仅供参考. //r_cname保存的是list类型,存储的是中文字符串 print r_cname //输出unicode编码格式,格式如下 [u'\u6e56\u4eba', u'\u7070\u718a', u'\u9ec4\u8702', u'\u70ed\u706

  • 用Python在Excel里画出蒙娜丽莎的方法示例

    之前看到过很多头条,说哪国某人坚持了多少年自学使用excel画画,效果十分惊艳. 对于他们的耐心我十分敬佩. 但是作为一个程序员,自然也得挑战一下自己. 这种需求,我们十分钟就可以完成! 基本思路 实现这个需求的基本思路是读取这张图片每一个像素的色彩值,然后给excel里的每一个单元格填充上颜色.所以主要用到的是PIL.openpyxl这两个库. PIL使用 PIL是Python里面做图像处理的时候十分常用的一个库,功能也是十分的强大,这里只需要用到PIL里一小部分的功能. from PIL i

随机推荐