In [94]:
import sys
sys.path.append("../codes/")
from Readfiles import getFnames
from DCdata import readReservoirDC_all
%pylab inline
In [95]:
directory = "../data/ChungCheonDC/"
In [96]:
fnames = getFnames(directory, dtype="apr", minimumsize=7000.)
In [97]:
fname_temp = fnames[0]
dat_temp,height_temp, ID = readReservoirDC_all(directory+fname_temp)
In [98]:
ntimes = len(fnames)
DATA = np.zeros((dat_temp.shape[0], ntimes))*np.nan
height = np.ones(ntimes)*np.nan
In [99]:
DATA.shape
Out[99]:
In [100]:
for i, fname in enumerate(fnames):
dat_temp,height_temp, ID = readReservoirDC_all(directory+fname)
if dat_temp.shape[0] == 380:
DATA[:,i] = dat_temp[:,-1]
height[i] = height_temp[0]
else:
print fname
In [101]:
a = ['1', '2', '3']
In [102]:
def strtofloat(input):
temp = ""
for i in input:
temp += i
return float(temp)
In [103]:
# dat_temp,height_temp, datalist = readReservoirDC_all(fnames[79])
# print fnames[79]
# # datalist = readReservoirDC_all(fnames[79])
In [132]:
print fnames[0]
In [105]:
locs = dat_temp[:,:4]
In [106]:
mida = locs[:,:2].sum(axis=1)
midb = locs[:,2:].sum(axis=1)
mid = (mida + midb)*0.5
dz = mida-midb
In [107]:
from ipywidgets import interact, IntSlider
In [108]:
from scipy import interpolate
In [109]:
x = np.linspace(mid.min(), mid.max(), 100)
z = np.linspace(dz.min(), dz.max(), 100)
# grid_x, grid_z = np.mgrid[np.min(mid):np.max(mid), np.min(dz):np.max(dz)]
grid_x, grid_z = np.meshgrid(x,z)
def vizDCtimeSeries(idatum, itime):
# idatum = 0
figsize(8,6)
fig = plt.figure()
ax1 = plt.subplot(211)
ax2 = plt.subplot(212)
# ax1.plot(mid, dz, '.')
grid_rho = griddata(mid, dz, DATA[:,itime], grid_x, grid_z, interp='linear')
grid_rho = grid_rho.reshape(grid_x.shape)
vmin, vmax = 50, 200.
ax1.contourf(grid_x, grid_z, grid_rho, 100, vmin =vmin, vmax = vmax, clim=(vmin, vmax), cmap="jet")
ax1.contourf(grid_x, grid_z, grid_rho, 100, vmin =vmin, vmax = vmax, clim=(vmin, vmax), cmap="jet")
ax1.scatter(mid, dz, s=20, c = DATA[:,itime], edgecolor="None", vmin =vmin, vmax = vmax, clim=(vmin, vmax))
# ax1.plot(grid_x.flatten(), grid_z.flatten(), 'k.')
ax1.plot(mid[idatum], dz[idatum], 'ro')
ax2.plot(DATA[idatum,:], 'k-', lw=2)
ax2.set_yscale('log')
ax2.set_ylim(vmin, vmax)
ax2_1 = ax2.twinx()
ax2_1.plot(height)
ax2_1.set_ylim(15, 21.)
ax2_1.plot(np.r_[itime, itime], np.r_[15, 21.], 'k--', lw=1)
ax1.text(0,0, fnames[itime])
In [ ]:
In [110]:
interact(vizDCtimeSeries, idatum=IntSlider(min=0, max=300, step=10, value=0),
itime=IntSlider(min=0, max=DATA.shape[1]-1, step=100, value=705))
Out[110]:
In [156]:
def vizDCtimeSeriesVariation(idatum, itime):
# idatum = 0
figsize(8,6)
fig = plt.figure()
ax1 = plt.subplot(211)
ax2 = plt.subplot(212)
# ax1.plot(mid, dz, '.')
itime_ref = 0
DATA_ref = DATA[:,itime_ref]
grid_rho = griddata(mid, dz, DATA[:,itime]/DATA_ref, grid_x, grid_z, interp='linear')
# grid_rho = griddata(mid, dz, (DATA[:,itime]-DATA_ref)/DATA_ref, grid_x, grid_z, interp='linear')
grid_rho = grid_rho.reshape(grid_x.shape)
vmin, vmax = 0.9, 1.1
ax1.contourf(grid_x, grid_z, grid_rho, 100, vmin =vmin, vmax = vmax, clim=(vmin, vmax), cmap="jet")
ax1.contourf(grid_x, grid_z, grid_rho, 100, vmin =vmin, vmax = vmax, clim=(vmin, vmax), cmap="jet")
ax1.scatter(mid, dz, s=20, c = DATA[:,itime]/DATA_ref, edgecolor="None", vmin =vmin, vmax = vmax, clim=(vmin, vmax))
# ax1.plot(grid_x.flatten(), grid_z.flatten(), 'k.')
ax1.plot(mid[idatum], dz[idatum], 'ro')
ax2.plot(DATA[idatum,:], 'k-', lw=2)
ax2.set_yscale('log')
vmin, vmax = 50., 200.
ax2.set_ylim(vmin, vmax)
ax2_1 = ax2.twinx()
ax2_1.plot(height)
ax2_1.set_ylim(15, 25.)
ax2_1.plot(np.r_[itime, itime], np.r_[15, 21.], 'k--', lw=1)
ax1.text(0,0, fnames[itime])
In [157]:
interact(vizDCtimeSeriesVariation, idatum=IntSlider(min=0, max=300, step=1, value=0),
itime=IntSlider(min=0, max=DATA.shape[1]-1, step=1, value=255))
Out[157]:
In [113]:
x = np.linspace(mid.min(), mid.max(), 100)
z = np.linspace(dz.min(), dz.max(), 40)
In [114]:
print z.min(), z.max()
print x.min(), x.max()
In [115]:
np.diff(x)
Out[115]:
In [116]:
from SimPEG import Mesh
from scipy import interpolate
In [117]:
hx = np.ones(110)*1.
hz = np.ones(40)*0.5
mesh2D = Mesh.TensorMesh([hx,hz], x0 = '0N')
In [118]:
len(range(90,900,4))
Out[118]:
In [119]:
timeind = range(0,1200,16)
hy = np.ones(len(timeind))
mesh = Mesh.TensorMesh([hx,hy,hz], x0 = '0N0')
itime_ref = 705
DATA_ref = DATA[:,itime_ref]
print "reference data", fnames[itime_ref]
# model = np.zeros((mesh2D.nC,len(timeind)))
# model_ratio = model.copy()
MODEL = np.zeros((mesh.nCx, mesh.nCy, mesh.nCz))
MODEL_ratio = np.zeros((mesh.nCx, mesh.nCy, mesh.nCz))
for i, itime in enumerate(timeind) :
print itime
F = interpolate.LinearNDInterpolator(np.c_[mid, dz], DATA[:,itime])
F_ratio = interpolate.LinearNDInterpolator(np.c_[mid, dz], abs(DATA[:,itime]-DATA_ref)/abs(DATA_ref))
MODEL[:,i,:] = F(mesh2D.gridCC).reshape((mesh.nCx, mesh.nCz), order="F")
MODEL_ratio[:,i,:] = F_ratio(mesh2D.gridCC).reshape((mesh.nCx, mesh.nCz), order="F")
In [120]:
from SimPEG import Utils
In [121]:
Model = Utils.mkvc(MODEL[:,:,:])
Model_ratio = Utils.mkvc(MODEL_ratio[:,:,:])
In [122]:
Model[np.isnan(Model)] = 1e-8
Model_ratio[np.isnan(Model_ratio)] = 1e-8
In [123]:
# mesh.writeModelUBC("/Users/sgkang/Dropbox/dammodel.txt", Model)
# mesh.writeModelUBC("/Users/sgkang/Dropbox/dammodel_ratio.txt", Model_ratio)
# mesh.writeUBC("/Users/sgkang/Dropbox/dammodel_mesh.txt")
mesh.writeModelUBC("dammodel705.txt", Model)
mesh.writeModelUBC("dammodel_ratio705.txt", Model_ratio)
mesh.writeUBC("dammodel_mesh705.txt")
In [124]:
plot(Model)
Out[124]:
In [125]:
mesh.plotSlice(Model, clim=(50., 200.))
Out[125]:
In [ ]:
In [126]:
dzu = np.unique(dz)
In [127]:
def profile_time(i_n, color):
figsize(6,3)
ind = np.argwhere(dz == dzu[::-1][i_n])
nskip = 5
for i in range(0,ind.size,nskip):
plt.semilogy(DATA[ind.flatten()[i],:], color)
plt.tight_layout()
plt.ylim(50, 200)
In [128]:
colors = ['k', 'b', 'r', 'g']
inds = [0, 2, 4, 6]
for i, ind in enumerate(inds):
profile_time(ind,colors[i])
In [129]:
profile_time(0, 'k')
In [130]:
profile_time(2, 'k')
In [131]:
profile_time(3, 'k')
In [ ]:
In [ ]: