Kriging插值计算

 

参考论文:      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]')