参考论文: http://people.ku.edu/~gbohling/cpe940
# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# Kriging.py
# Created on: 2014-06-12 10:14:21.00000
# (generated by ArcGIS/ModelBuilder)
# Description:
# --------------------------------------------------------------------------- # Import arcpy module
import os
import math
import sys from pylab import *
import numpy as np
from pandas import DataFrame, Series
from scipy.spatial.distance import pdist, squareform # 计算距离 # dataXYV = [{"x":12100,"y":8300,"v":14.6515},{"x":5300,"y":8700,"v":14.5093},{"x":3500,"y":13900,"v":14.0639},{"x":5100,"y":1900,"v":15.1084},{"x":9900,"y":13700,"v":13.919},{"x":2900,"y":900,"v":13.1304},{"x":7900,"y":6700,"v":14.5724},{"x":16900,"y":4900,"v":15.0814},{"x":18700,"y":1500,"v":13.91},{"x":2700,"y":2100,"v":13.4024},{"x":10700,"y":5100,"v":14.9395},{"x":7500,"y":12900,"v":15.2159},{"x":5500,"y":11100,"v":14.5777},{"x":9500,"y":9100,"v":14.2483},{"x":15300,"y":3100,"v":14.4281},{"x":4700,"y":9700,"v":15.2606},{"x":16700,"y":15700,"v":16.1859},{"x":19500,"y":9700,"v":14.2079},{"x":16900,"y":13100,"v":16.9583},{"x":900,"y":3700,"v":13.8354},{"x":500,"y":11900,"v":14.1859},{"x":9100,"y":1300,"v":14.0381},{"x":9100,"y":13700,"v":14.3685},{"x":9900,"y":12900,"v":13.4018},{"x":6300,"y":100,"v":15.8953},{"x":3700,"y":5100,"v":12.8667},{"x":16300,"y":900,"v":15.1039},{"x":18300,"y":13500,"v":15.7736},{"x":9500,"y":6900,"v":14.1333},{"x":17900,"y":3100,"v":13.3369},{"x":9900,"y":15500,"v":15.1362},{"x":7100,"y":8900,"v":15.0847},{"x":19300,"y":7100,"v":14.2498},{"x":2300,"y":5700,"v":12.6811},{"x":7300,"y":8900,"v":14.9384},{"x":13900,"y":3700,"v":15.6005},{"x":8500,"y":10100,"v":13.7796},{"x":8100,"y":8700,"v":15.2907},{"x":14700,"y":11900,"v":15.6881},{"x":6300,"y":2300,"v":15.3677},{"x":11900,"y":12900,"v":14.3283},{"x":18100,"y":7100,"v":14.7374},{"x":11300,"y":7100,"v":15.0547},{"x":12500,"y":3100,"v":14.8889},{"x":2700,"y":12700,"v":14.436},{"x":2700,"y":4300,"v":12.1491},{"x":8500,"y":11300,"v":13.624},{"x":1500,"y":900,"v":14.188},{"x":7300,"y":1300,"v":14.9072},{"x":10700,"y":4100,"v":15.2029},{"x":7100,"y":1900,"v":15.3468},{"x":3900,"y":8500,"v":15.939},{"x":17100,"y":6100,"v":15.7269},{"x":14100,"y":10100,"v":15.3238},{"x":11500,"y":4900,"v":14.0445},{"x":13300,"y":15700,"v":14.4032},{"x":1900,"y":12100,"v":14.3586},{"x":15100,"y":2900,"v":14.6007},{"x":6500,"y":900,"v":16.1458},{"x":8900,"y":6100,"v":15.7727},{"x":4500,"y":2300,"v":13.6234},{"x":12900,"y":10300,"v":15.1024},{"x":10900,"y":5700,"v":15.3546},{"x":3500,"y":700,"v":13.8431},{"x":16300,"y":3700,"v":14.9427},{"x":900,"y":5100,"v":14.4139},{"x":12900,"y":12900,"v":13.6177},{"x":15300,"y":9300,"v":16.3787},{"x":7300,"y":6900,"v":14.258},{"x":16300,"y":12500,"v":15.7772},{"x":100,"y":8900,"v":14.6553},{"x":1700,"y":11700,"v":14.3627},{"x":17500,"y":11100,"v":15.9659},{"x":14900,"y":8300,"v":16.0095},{"x":8300,"y":10900,"v":13.9639},{"x":4100,"y":14500,"v":14.2649},{"x":11100,"y":15300,"v":15.7684},{"x":500,"y":4900,"v":14.591},{"x":13100,"y":1500,"v":15.1377},{"x":18900,"y":1700,"v":14.095},{"x":3500,"y":7500,"v":15.1486},{"x":3700,"y":6900,"v":13.9584},{"x":14500,"y":13300,"v":14.7381},{"x":4900,"y":9100,"v":15.0689},{"x":9700,"y":5700,"v":15.8042}]
dataXYV = [{"x":12100.00,"y":8300.00,"v":14.6515},
{"x":5300.00,"y":8700.00,"v":14.5093},
{"x":3500.00,"y":13900.00,"v":14.0639},
{"x":5100.00,"y":1900.00,"v":15.1084},
{"x":9900.00,"y":13700.00,"v":13.919},
{"x":2900.00,"y":900.00,"v":13.1304},
{"x":7900.00,"y":6700.00,"v":14.5724},
{"x":16900.00,"y":4900.00,"v":15.0814},
{"x":18700.00,"y":1500.00,"v":13.91},
{"x":2700.00,"y":2100.00,"v":13.4024},
{"x":10700.00,"y":5100.00,"v":14.9395},
{"x":7500.00,"y":12900.00,"v":15.2159},
{"x":5500.00,"y":11100.00,"v":14.5777},
{"x":9500.00,"y":9100.00,"v":14.2483},
{"x":15300.00,"y":3100.00,"v":14.4281},
{"x":4700.00,"y":9700.00,"v":15.2606},
{"x":16700.00,"y":15700.00,"v":16.1859},
{"x":19500.00,"y":9700.00,"v":14.2079},
{"x":16900.00,"y":13100.00,"v":16.9583},
{"x":900.00,"y":3700.00,"v":13.8354},
{"x":500.00,"y":11900.00,"v":14.1859},
{"x":9100.00,"y":1300.00,"v":14.0381},
{"x":9100.00,"y":13700.00,"v":14.3685},
{"x":9900.00,"y":12900.00,"v":13.4018},
{"x":6300.00,"y":100.00,"v":15.8953},
{"x":3700.00,"y":5100.00,"v":12.8667},
{"x":16300.00,"y":900.00,"v":15.1039},
{"x":18300.00,"y":13500.00,"v":15.7736},
{"x":9500.00,"y":6900.00,"v":14.1333},
{"x":17900.00,"y":3100.00,"v":13.3369},
{"x":9900.00,"y":15500.00,"v":15.1362},
{"x":7100.00,"y":8900.00,"v":15.0847},
{"x":19300.00,"y":7100.00,"v":14.2498},
{"x":2300.00,"y":5700.00,"v":12.6811},
{"x":7300.00,"y":8900.00,"v":14.9384},
{"x":13900.00,"y":3700.00,"v":15.6005},
{"x":8500.00,"y":10100.00,"v":13.7796},
{"x":8100.00,"y":8700.00,"v":15.2907},
{"x":14700.00,"y":11900.00,"v":15.6881},
{"x":6300.00,"y":2300.00,"v":15.3677},
{"x":11900.00,"y":12900.00,"v":14.3283},
{"x":18100.00,"y":7100.00,"v":14.7374},
{"x":11300.00,"y":7100.00,"v":15.0547},
{"x":12500.00,"y":3100.00,"v":14.8889},
{"x":2700.00,"y":12700.00,"v":14.436},
{"x":2700.00,"y":4300.00,"v":12.1491},
{"x":8500.00,"y":11300.00,"v":13.624},
{"x":1500.00,"y":900.00,"v":14.188},
{"x":7300.00,"y":1300.00,"v":14.9072},
{"x":10700.00,"y":4100.00,"v":15.2029},
{"x":7100.00,"y":1900.00,"v":15.3468},
{"x":3900.00,"y":8500.00,"v":15.939},
{"x":17100.00,"y":6100.00,"v":15.7269},
{"x":14100.00,"y":10100.00,"v":15.3238},
{"x":11500.00,"y":4900.00,"v":14.0445},
{"x":13300.00,"y":15700.00,"v":14.4032},
{"x":1900.00,"y":12100.00,"v":14.3586},
{"x":15100.00,"y":2900.00,"v":14.6007},
{"x":6500.00,"y":900.00,"v":16.1458},
{"x":8900.00,"y":6100.00,"v":15.7727},
{"x":4500.00,"y":2300.00,"v":13.6234},
{"x":12900.00,"y":10300.00,"v":15.1024},
{"x":10900.00,"y":5700.00,"v":15.3546},
{"x":3500.00,"y":700.00,"v":13.8431},
{"x":16300.00,"y":3700.00,"v":14.9427},
{"x":900.00,"y":5100.00,"v":14.4139},
{"x":12900.00,"y":12900.00,"v":13.6177},
{"x":15300.00,"y":9300.00,"v":16.3787},
{"x":7300.00,"y":6900.00,"v":14.258},
{"x":16300.00,"y":12500.00,"v":15.7772},
{"x":100.00,"y":8900.00,"v":14.6553},
{"x":1700.00,"y":11700.00,"v":14.3627},
{"x":17500.00,"y":11100.00,"v":15.9659},
{"x":14900.00,"y":8300.00,"v":16.0095},
{"x":8300.00,"y":10900.00,"v":13.9639},
{"x":4100.00,"y":14500.00,"v":14.2649},
{"x":11100.00,"y":15300.00,"v":15.7684},
{"x":500.00,"y":4900.00,"v":14.591},
{"x":13100.00,"y":1500.00,"v":15.1377},
{"x":18900.00,"y":1700.00,"v":14.095},
{"x":3500.00,"y":7500.00,"v":15.1486},
{"x":3700.00,"y":6900.00,"v":13.9584},
{"x":14500.00,"y":13300.00,"v":14.7381},
{"x":4900.00,"y":9100.00,"v":15.0689},
{"x":9700.00,"y":5700.00,"v":15.8042}] length = len(dataXYV)
distanceMatrix = [[] for i in range(length)]
index = 0
distTotal = 0;
distMin = 1.0e15
distMax = -1.0e15
distAver = 0
for x in dataXYV:
for y in dataXYV:
z = math.sqrt((x['x']-y['x'])*(x['x']-y['x'])+(x['y']-y['y'])*(x['y']-y['y']))
distTotal += z
if z > distMax:
distMax = z
if z < distMin and x != y:
distMin = z
distanceMatrix[index].append(z)
index += 1
distAver = distTotal / (length * length - length)
dataInfo = {'count':(length * length - length), 'distAver':distAver, 'distMin':distMin, 'distMax':distMax}
#print dataInfo '''
for i in range(0, length):
for j in range(0, length):
print(int(lists[i][j])),
print(';')
''' '''
查找点对
'''
def findPairs(dataXYV, distanceMatrix, minValue, maxValue):
totalDistance = 0;
count = 0;
minDistance = 1.0e15
maxDistance = -1.0e15
averageDistance = 0
pairs = []
for i in range(0, length):
for j in range(i+1, length):
if distanceMatrix[i][j]>minValue and distanceMatrix[i][j]<=maxValue:
#if math.fabs(dataXYV[i]['x']-dataXYV[j]['x'])>minValue and math.fabs(dataXYV[i]['y']-dataXYV[j]['y'])>minValue and math.fabs(dataXYV[i]['x']-dataXYV[j]['x'])<=maxValue and math.fabs(dataXYV[i]['y']-dataXYV[j]['y'])<=maxValue:
# print(int(lists[i][j])),
totalDistance += distanceMatrix[i][j]
count += 1
if distanceMatrix[i][j] >= maxDistance:
maxDistance = distanceMatrix[i][j]
if distanceMatrix[i][j] <= minDistance:
minDistance = distanceMatrix[i][j]
pair = {'i':i,'j':j,'iv':dataXYV[i]['v'],'jv':dataXYV[j]['v'],'dist':distanceMatrix[i][j]}
pairs.append(pair)
#print(count)
averageDistance = totalDistance / count
info = {'count':count, 'distAver':averageDistance, 'distMin':minDistance, 'distMax':maxDistance}
#print info
#print pairs
return pairs '''
计算统计信息: 协方差、相关系数、半方变异
'''
def computeStatisticInfo(pairs):
pairCount = len(pairs)
distanceTotal = 0
distanceAverage = 0
#
v1v2Total=0
v1Total=0
v2Total=0
#
v1v1Total=0
v2v2Total=0
#
v1_v2Total=0
#
for x in pairs:
val1 = x['iv']
val2 = x['jv']
distanceTotal = distanceTotal + x['dist']
v1v2Total = v1v2Total + val1 * val2
v1Total = v1Total + val1
v2Total = v2Total + val2
#
v1v1Total = v1v1Total + val1 * val1
v2v2Total = v2v2Total + val2 * val2
#
v1_v2Total = v1_v2Total + math.pow(val1 - val2, 2)
#
distanceAverage = distanceTotal / pairCount
v1v2Covariance = v1v2Total / pairCount - v1Total * v2Total / (pairCount * pairCount)
v1v2Corelation = (v1v2Total*pairCount - v1Total * v2Total) / math.sqrt(v1v1Total * pairCount - v1Total * v1Total) / math.sqrt(v2v2Total * pairCount - v2Total * v2Total)
v1v2Semivariance = v1_v2Total / (pairCount * 2)
statisticInfo = {'covariance':v1v2Covariance, 'corelation':v1v2Corelation, 'semivariance':v1v2Semivariance, 'count':pairCount, 'distAver':distanceAverage} #
# print statisticInfo
return statisticInfo '''
计算各种lagSize下的统计信息: 协方差、相关系数、半方变异
'''
def staticInfoAll(dataXYV, distanceMatrix, lagCellSize, lagCount):
semiLagCellSize = lagCellSize / 2
pairsStaticInfos = []
for i in range(0, lagCount-1): lagSize = lagCellSize * i
lagMin = lagSize - semiLagCellSize
if lagMin < 0:
lagMin = 0
lagMax = lagSize + semiLagCellSize
#print(lagMin, lagMax, lagSize) '''
lagMin = lagCellSize * i
lagMax = lagCellSize * (i + 1)
lagSize = (lagCellSize * i + lagCellSize * (i + 1))/2
#print lagMin, lagMax, lagSize
'''
pairs = findPairs(dataXYV, distanceMatrix, lagMin, lagMax)
statisticInfo = computeStatisticInfo(pairs)
statisticInfo['lagSize'] = lagSize
print(lagMin, lagMax, statisticInfo['lagSize'], statisticInfo['count'], statisticInfo['distAver'], statisticInfo['covariance'], statisticInfo['corelation'], statisticInfo['semivariance'])
pairsStaticInfos.append(statisticInfo)
return pairsStaticInfos
'''
def computeC0(dataXYV):
valueTotal = 0;
valueAver = 0;
variance = 0;
for x in dataXYV:
valueTotal += x['v']
valueAver = valueTotal / length
print valueAver
for x in dataXYV:
variance += math.pow((x['v'] - valueAver),2)
variance = variance / length
# variance = math.sqrt(variance)
print variance computeC0(dataXYV)
''' def optimization(pairsStaticInfos):
aMin = 0
cMin = 0
varianceTotalMin = 1.0e45
hSize = len(pairsStaticInfos)
for aValue in range( 3500, 4500):
for c1Value in range(60, 99):
cValue = c1Value / 100.0
#cValue = 0.78
#print(aValue, cValue)
#print(aValue, cValue)
varianceTotal = 0
for x in pairsStaticInfos:
y = spherical( x['distAver'], aValue, cValue ) # distAver lagSize
#print y, x['semivariance']
varianceTotal = varianceTotal + ((y - x['semivariance'])**2.0) * x['count']
varianceTotal = varianceTotal / hSize
#print varianceTotal
if varianceTotal <= varianceTotalMin:
varianceTotalMin = varianceTotal
aMin = aValue
cMin = cValue
#print(aMin, cMin, varianceTotalMin) para = {"a":aMin, "c0":cMin}
print(para)
return para def spherical( h, a, C0):
if h <= a:
return (C0*( 1.5*h/a - 0.5*(h/a)**3.0 ))
else:
return C0 pairsStaticInfos = staticInfoAll(dataXYV, distanceMatrix, 1000, 12) # 500, 25 para = optimization(pairsStaticInfos) lagSize = []
semivariance = []
modelY = []
modelYY = []
#
lagSize.append(0)
semivariance.append(0)
modelY.append(0)
for x in pairsStaticInfos:
y = spherical( x['distAver'], para['a'], para['c0'] )
#yy= spherical( x['distAver'], 4141, para['c0'] )
lagSize.append(x['distAver'])
semivariance.append(x['semivariance'])
modelY.append(y)
modelYY.append(yy) plot( lagSize, semivariance, 'o' )
plot( lagSize, modelY, '.-' ) ;
title('Spherical Model')
ylabel('Semivariance')
xlabel('Lag [m]')