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]')
