def _many_sample(self, xloc, yloc, xa, ya, vra):
"Solve the Kriging System with more than one sample"
na = len(vra)
# number of equations, for simple kriging there're na,
# for ordinary there're na + 1
neq = na + self.ktype
# Establish left hand side covariance matrix:
left = np.full((neq, neq), np.nan)
for i, j in product(xrange(na), xrange(na)):
if np.isnan(left[j, i]):
left[i, j] = self._cova2(xa[i], ya[i], xa[j], ya[j])
else:
left[i, j] = left[j, i]
# Establish the Right Hand Side Covariance:
right = list()
for j in xrange(na):
xx = xa[j] - xloc
yy = ya[j] - yloc
if not self.block_kriging:
cb = self._cova2(xx, yy, self.xdb[0], self.ydb[0])
else:
cb = 0.0
for i in xrange(self.ndb):
cb += self._cova2(xx, yy, self.xdb[i], self.ydb[i])
dx = xx - self.xdb[i]
dy = yy - self.ydb[i]
if dx*dx + dy*dy < np.finfo('float').eps:
cb -= self.c0
cb /= self.ndb
right.append(cb)
if self.ktype == 1: # for ordinary kriging
# Set the unbiasedness constraint
left[neq-1, :-1] = self.unbias
left[:-1, neq-1] = self.unbias
left[-1, -1] = 0
right.append(self.unbias)
# Solve the kriging system
s = None
try:
s = linalg.solve(left, right)
except linalg.LinAlgError as inst:
print("Warning kb2d: singular matrix for block " + \
"{},{}".format((xloc-self.xmn)/self.xsiz,
(yloc-self.ymn)/self.ysiz))
return np.nan, np.nan
estv = self.block_covariance
if self.ktype == 1: # ordinary kriging
estv -= s[-1]*self.unbias # s[-1] is mu
est = np.sum(s[:na]*vra[:na])
estv -= np.sum(s[:na]*right[:na])
if self.ktype == 0: # simple kriging
est += (1 - np.sum(s[:na])) * self.skmean
return est, estv
评论列表
文章目录
