In [9]:
import sys
sys.path.append("../codes/")
from Readfiles import getFnames
from DCdata import readReservoirDC_all
%pylab inline
In [10]:
directory = "../data/ChungCheonDC/20151218000000.apr"
In [11]:
fnames = getFnames(directory, dtype="apr", minimumsize=7000.)
In [12]:
fname_temp = fnames[0]
dat_temp,height_temp = readReservoirDC_all(fname_temp)
In [ ]:
ntimes = len(fnames)
DATA = np.zeros((dat_temp.shape[0], ntimes))*np.nan
height = np.ones(ntimes)*np.nan
In [ ]:
DATA.shape
In [ ]:
for i, fname in enumerate(fnames):
dat_temp,height_temp = readReservoirDC_all(fname)
if dat_temp.shape[0] == 380:
DATA[:,i] = dat_temp[:,-1]
height[i] = height_temp[0]
else:
print fname
In [ ]:
a = ['1', '2', '3']
In [ ]:
def strtofloat(input):
temp = ""
for i in input:
temp += i
return float(temp)
In [ ]:
# dat_temp,height_temp, datalist = readReservoirDC_all(fnames[79])
# print fnames[79]
# # datalist = readReservoirDC_all(fnames[79])
In [ ]:
print fnames[79]
In [ ]:
locs = dat_temp[:,:4]
In [ ]:
mida = locs[:,:2].sum(axis=1)
midb = locs[:,2:].sum(axis=1)
mid = (mida + midb)*0.5
dz = mida-midb
In [ ]:
from ipywidgets import interact, IntSlider
In [ ]:
from scipy import interpolate
In [ ]:
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 [ ]:
interact(vizDCtimeSeries, idatum=IntSlider(min=0, max=300, step=10, value=0),
itime=IntSlider(min=0, max=DATA.shape[1]-1, step=100, value=0))
In [ ]:
In [ ]:
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 = 790
DATA_ref = DATA[:,itime_ref]
grid_rho = griddata(mid, dz, DATA[:,itime]/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, 21.)
ax2_1.plot(np.r_[itime, itime], np.r_[15, 21.], 'k--', lw=1)
ax1.text(0,0, fnames[itime])
In [ ]:
interact(vizDCtimeSeriesVariation, idatum=IntSlider(min=0, max=300, step=10, value=0),
itime=IntSlider(min=0, max=DATA.shape[1]-1, step=4, value=0))
In [ ]:
DATA.shape
In [ ]:
x = np.linspace(mid.min(), mid.max(), 100)
z = np.linspace(dz.min(), dz.max(), 40)
In [ ]:
print z.min(), z.max()
print x.min(), x.max()
In [ ]:
np.diff(x)
In [ ]:
from SimPEG import Mesh
from scipy import interpolate
In [ ]:
hx = np.ones(110)*1.
hy = np.ones(40)*0.5
mesh2D = Mesh.TensorMesh([hx,hy], x0 = '0N')
In [ ]:
len(range(90,900,4))
In [ ]:
timeind = range(0,1200,16)
model = np.zeros((mesh2D.nC,len(timeind)))
hz = np.ones(len(timeind))
mesh = Mesh.TensorMesh([hx,hy,hz], x0 = '0N0')
itime_ref = 790
DATA_ref = DATA[:,itime_ref]
model_ratio = model.copy()
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_ratio[:,i] = F_ratio(mesh2D.gridCC)
model[:,i] = F(mesh2D.gridCC)
In [ ]:
model.shape
In [ ]:
from SimPEG import Utils
In [ ]:
Model = Utils.mkvc(model[:,::-1])
Model_ratio = Utils.mkvc(model_ratio[:,::-1])
In [ ]:
Model[np.isnan(Model)] = 1e-8
Model_ratio[np.isnan(Model_ratio)] = 1e-8
In [ ]:
mesh.writeModelUBC("./dammodel.txt", Model)
mesh.writeModelUBC("./dammodel_ratio.txt", Model_ratio)
In [ ]:
mesh.writeUBC("./dammodel_mesh.txt")
In [ ]:
!dir
In [ ]:
plot(Model)
In [ ]:
mesh.plotSlice(Model, clim=(50., 200.))
In [ ]:
In [ ]:
dzu = np.unique(dz)
In [ ]:
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 [ ]:
colors = ['k', 'b', 'r', 'g']
inds = [0, 2, 4, 6]
for i, ind in enumerate(inds):
profile_time(ind,colors[i])
In [ ]:
profile_time(0, 'k')
In [ ]:
profile_time(2, 'k')
In [ ]:
profile_time(3, 'k')
In [ ]:
profile_time(4, 'k')
In [ ]:
profile_time(5, 'k')
In [ ]:
profile_time(6, 'k')
In [ ]:
profile_time(7, 'k')
In [ ]: