In [1]:
import sys
sys.path.append("../codes/")
from Readfiles import getFnames
from DCdata import readReservoirDC_all
%pylab inline
In [2]:
directory = "../data/ChungCheonDC/"
In [3]:
fnames = getFnames(directory, dtype="apr", minimumsize=7000.)
In [4]:
fname_temp = fnames[0]
dat_temp,height_temp, ID = readReservoirDC_all(directory+fname_temp)
In [5]:
ntimes = len(fnames)
DATA = np.zeros((dat_temp.shape[0], ntimes))*np.nan
height = np.ones(ntimes)*np.nan
In [6]:
len(fnames)
Out[6]:
In [7]:
fnamesactive = []
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]
# print fname
fnamesactive.append(fname)
else:
print fname
In [8]:
a = ['1', '2', '3']
In [9]:
def strtofloat(input):
temp = ""
for i in input:
temp += i
return float(temp)
In [10]:
# dat_temp,height_temp, datalist = readReservoirDC_all(fnames[79])
# print fnames[79]
# # datalist = readReservoirDC_all(fnames[79])
In [11]:
print fnames[79]
In [12]:
locs = dat_temp[:,:4]
In [13]:
mida = locs[:,:2].sum(axis=1)
midb = locs[:,2:].sum(axis=1)
mid = (mida + midb)*0.5
dz = mida-midb
In [14]:
from ipywidgets import interact, IntSlider
In [15]:
from scipy import interpolate
In [16]:
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 [17]:
interact(vizDCtimeSeries, idatum=IntSlider(min=0, max=300, step=5, value=0),
itime=IntSlider(min=0, max=DATA.shape[1]-1, step=100, value=0))
Out[17]:
In [ ]:
In [18]:
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 [19]:
interact(vizDCtimeSeriesVariation, idatum=IntSlider(min=0, max=300, step=1, value=0),
itime=IntSlider(min=0, max=DATA.shape[1]-1, step=20, value=0))
Out[19]:
In [20]:
DATA.shape
Out[20]:
In [21]:
x = np.linspace(mid.min(), mid.max(), 100)
z = np.linspace(dz.min(), dz.max(), 40)
In [22]:
print z.min(), z.max()
print x.min(), x.max()
In [23]:
np.diff(x)
Out[23]:
In [24]:
from SimPEG import Mesh
from scipy import interpolate
In [25]:
hx = np.ones(110)*1.
hz = np.ones(40)*0.5
mesh2D = Mesh.TensorMesh([hx,hz], x0 = '0N')
In [28]:
len(range(90,900,4))
Out[28]:
In [27]:
fnamesfinal = []
for ind in timeind:
fnamesfinal.append(fnamesactive[ind])
In [ ]:
timeind = range(0,1200,4)
hy = np.ones(len(timeind))
mesh = Mesh.TensorMesh([hx,hy,hz], x0 = '00N')
itime_ref = 790
DATA_ref = DATA[:,itime_ref]
print "reference data", fnamesactive[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 [ ]:
MODEL.shape
In [ ]:
from SimPEG import Utils
In [ ]:
Model = Utils.mkvc(MODEL[:,:,:])
Model_ratio = Utils.mkvc(MODEL_ratio[:,:,:])
In [ ]:
Model[np.isnan(Model)] = 1e-8
Model_ratio[np.isnan(Model_ratio)] = 1e-8
In [ ]:
# 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("dammodel901.txt", Model)
mesh.writeModelUBC("dammodel901_ratio.txt", Model_ratio)
mesh.writeUBC("dammodel_mesh901.txt")
In [ ]:
iloc = 10
fnamesfinal[iloc]
In [ ]:
plot(Model)
In [ ]:
In [ ]:
mesh.plotSlice(Model, clim=(50., 200.))
In [ ]:
In [74]:
dzu = np.unique(dz)
In [75]:
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 [76]:
colors = ['k', 'b', 'r', 'g']
inds = [0, 2, 4, 6]
for i, ind in enumerate(inds):
profile_time(ind,colors[i])
In [77]:
profile_time(0, 'k')
In [78]:
profile_time(2, 'k')
In [79]:
profile_time(3, 'k')
In [80]:
profile_time(4, 'k')
In [81]:
profile_time(5, 'k')
In [82]:
profile_time(6, 'k')
In [83]:
profile_time(7, 'k')
In [ ]:
In [ ]: