Closed Ml0749276 closed 1 year ago
Hey there!
I got this while I tried to OrdinaryKriging in a loop.
Where did it say something about SVD as you have written in the title?
The error output is as below
Adjusting data for anisotropy... Initializing variogram model...
Is there really an error? Or is it hanging at this point? If it is hanging, I think it could be at this line: https://github.com/GeoStat-Framework/PyKrige/blob/b33e5d87fb26abf0c919018bddbed65156d2ed37/pykrige/core.py#L455
Or did you miss to copy the full output? How big is the size of x
and y
in you case?
Cheers, Sebastian
Hi Sebastian,
Here is the full script
# filecoding=utf-8
# 用于对测点数据进行插值,得到平面值
import math
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np
from pykrige.ok import OrdinaryKriging
import readZdatafile
import sitesclass
D2R = math.pi/180
R2D = 180/math.pi
def DegtoFirstQuadrant(angle):
result = math.atan(math.tan(angle*D2R))*R2D
if result < 0:
result = result+180
return result
Tsite = sitesclass.Tsite
Tsites = sitesclass.Tsites
ObseSites = Tsites()
obsesitelist = []
obsesitelist.append('0m0d5%.ori_dat')
obsesitelist.append('0m85d5%.ori_dat')
# 待绘制频率
freqlist = [0.25, 0.05, 0.025, 0.0125, 0.00625]
for i in range(len(obsesitelist)):
obsedir = obsesitelist[i]
readZdatafile.readfiles(obsedir, ObseSites)
# 设置网格
xstep = (ObseSites.xmax-ObseSites.xmin+20)/(len(ObseSites.items)-1)
ystep = (ObseSites.ymax-ObseSites.ymin+20)/(len(ObseSites.items)-1)
gridx = np.arange(ObseSites.xmin/1000.-10, ObseSites.xmax/1000.+10, min(xstep, ystep)/1000.)
gridy = np.arange(ObseSites.ymin/1000.-10, ObseSites.ymax/1000.+10, min(xstep, ystep)/1000.)
fig, ax = plt.subplots(4, 5, dpi=180, sharex='all', sharey='all', figsize=(4, 5))
# 循环频率,绘制观测数据
for index in range(len(freqlist)):
freq = freqlist[index]
x = np.zeros(len(ObseSites.items))
y = np.zeros(len(ObseSites.items))
rxy = np.zeros(len(ObseSites.items))
pxy = np.zeros(len(ObseSites.items))
ryx = np.zeros(len(ObseSites.items))
pyx = np.zeros(len(ObseSites.items))
component = [rxy, ryx, pxy, pyx]
for obse in ObseSites.items:
iobse = ObseSites.items.index(obse)
for j in range(len(obse.freqs)):
if math.isclose(obse.freqs[j], freq):
x[iobse] = obse.x/1000.
y[iobse] = obse.y/1000.
rxy[iobse] = math.log10(obse.rxy[j].re)
pxy[iobse] = DegtoFirstQuadrant(obse.pxy[j].re)
ryx[iobse] = math.log10(obse.ryx[j].re)
pyx[iobse] = DegtoFirstQuadrant(obse.pyx[j].re)
break
for j in range(len(component)):
if j < 2:
vmin = 1
vmax = 5
else:
vmin = 0
vmax = 100
ok = OrdinaryKriging(y, x, component[j], variogram_model='exponential', verbose=True)
z, ss = ok.execute('grid', gridy, gridx)
ax[j, index].pcolormesh(gridy, gridx, z, cmap='jet_r',
vmin=vmin, vmax=vmax, linewidth=0)
ax[j, index].scatter(y, x, c='black', s=5)
ax[j, index].set_aspect(1)
ax[j, index].tick_params(labelsize=5)
# 调整坐标轴
ax[j, 0].set_ylabel('X(km)', fontsize=5)
# 设置图头
ax[0, index].set_title(str(freq)+'Hz', fontsize=5, pad=0.6)
# 调整坐标轴
ax[3, index].set_xlabel('Y(km)', fontsize=5)
# 调整布局
fig.subplots_adjust(hspace=0.1, wspace=0.05, right=0.85,
left=0.11, top=0.95, bottom=0.1)
# 设置右侧分量名
fig.text(0.85, 0.85, 'Rxy', fontsize=5, va='center', rotation=90)
fig.text(0.85, 0.63, 'Ryx', fontsize=5, va='center', rotation=90)
fig.text(0.85, 0.41, 'Pxy', fontsize=5, va='center', rotation=90)
fig.text(0.85, 0.20, 'Pyx', fontsize=5, va='center', rotation=90)
# 设置色标
# 视电阻率
norm = mpl.colors.Normalize(vmin=1, vmax=5)
m = cm.ScalarMappable(norm=norm, cmap=cm.jet_r)
pos = fig.add_axes([0.9, 0.55, 0.02, 0.35]) # 位置[左、下、宽度、高度]
cb = plt.colorbar(m, cax=pos)
cb.ax.tick_params(labelsize=5)
cb.ax.set_title('log10(rho)', fontsize=5, loc='left')
# 相位
norm = mpl.colors.Normalize(vmin=0, vmax=100)
m = cm.ScalarMappable(norm=norm, cmap=cm.jet_r)
pos = fig.add_axes([0.9, 0.11, 0.02, 0.35]) # 位置[左、下、宽度、高度]
cb = plt.colorbar(m, cax=pos)
cb.ax.tick_params(labelsize=5)
cb.ax.set_title('degree', fontsize=5, loc='left')
plt.show()
plt.close()
I write this script to genarate figures, the size of 'x' and 'y' is 33, and the size of the mesh is [38,87]. During the whole workflow, I read two data files and generate two figures for each one, after the first figure generated successfully, the second one will jump out the error 'SVD did not converge' as attached picuture. In the beginning, I thought it might be the problem of my datafile, so I switched the order of the two datafile and read the 2nd first, then it worked and generated figure for the 2nd data file, but the same error occured while processing the 1st data file. So I got confused. Can you help me solving this? Thank you!
Attched are the error message captured and two figures generated. PyKrige is a convenient tool for the interpolation and plotting!
The error message
The 1st data file's figure
The 2nd data file's figure
And another question, this subfigure seems not interpolate correctly, do you have some ideas about this? I can send you all the data if you want, this situation not only appeared here, but many other places.
Also in you case it could help to use the pseudo-inverse version from this branch: https://github.com/GeoStat-Framework/PyKrige/tree/pseudo_inv
Then you can alter your call:
ok = OrdinaryKriging(y, x, component[j], variogram_model='exponential', verbose=True, pseudo_inv=True)
Maybe this helps. This will be in the next version.
Hi Sebastian,
Thanks for your reply, I have tried as you said, it makes no difference.
Could you check your data for nan
s and inf
s? That is often related to a svd problem
Hi Sebastian,
I am on a business trip these days, I check the data and no nan
or inf
is found.
Closing due to inactivity.
I got this while I tried to OrdinaryKriging in a loop.
Here is part of code
The error output is as below