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]:
DATA.shape
Out[6]:
In [7]:
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 [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[0]
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 [ ]:
In [17]:
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[17]:
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 = 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 [19]:
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[19]:
In [20]:
x = np.linspace(mid.min(), mid.max(), 100)
z = np.linspace(dz.min(), dz.max(), 40)
In [21]:
print z.min(), z.max()
print x.min(), x.max()
In [22]:
np.diff(x)
Out[22]:
In [23]:
from SimPEG import Mesh
from scipy import interpolate
In [24]:
hx = np.ones(110)*1.
hz = np.ones(40)*0.5
mesh2D = Mesh.TensorMesh([hx,hz], x0 = '0N')
In [25]:
len(range(90,900,4))
Out[25]:
In [26]:
timeind = range(0,1200,16)
hy = np.ones(len(timeind))
mesh = Mesh.TensorMesh([hx,hy,hz], x0 = '0N0')
itime_ref = 922 ## 705 : reference data 20150721000000.apr
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 [27]:
from SimPEG import Utils
In [28]:
Model = Utils.mkvc(MODEL[:,:,:])
Model_ratio = Utils.mkvc(MODEL_ratio[:,:,:])
In [29]:
Model[np.isnan(Model)] = 1e-8
Model_ratio[np.isnan(Model_ratio)] = 1e-8
In [30]:
# 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("dammodel913.txt", Model)
mesh.writeModelUBC("dammodel_ratio913.txt", Model_ratio)
mesh.writeUBC("dammodel_mesh913.txt")
In [31]:
plot(Model)
Out[31]:
In [32]:
mesh.plotSlice(Model, clim=(50., 200.))
Out[32]:
In [ ]:
In [33]:
dzu = np.unique(dz)
In [34]:
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 [35]:
colors = ['k', 'b', 'r', 'g']
inds = [0, 2, 4, 6]
for i, ind in enumerate(inds):
profile_time(ind,colors[i])
In [36]:
profile_time(0, 'k')
In [37]:
profile_time(6, 'k')
In [38]:
import numpy as np
import scipy.stats as st
In [39]:
n = 100 # length of the data
x = np.random.randn(n)
y = 3 + 7*x + np.random.randn(n)
#y = 3 + 7*x + np.random.normal(n)
#y = state_means
In [40]:
b, a, r, p, e = st.linregress(x, y)
In [41]:
eps = y - a - b*x # error of fitting and measured data
x1 = np.linspace(0, 1) # x axis to plot the PI
# variace of fitting error
e_pi = np.var(eps)*(1+1.0/n + (x1-x.mean())**2/np.sum((x-x.mean())**2))
# z value using the t distribution and with dof = n-2
z = st.t.ppf(0.95, n-2)
# prediction interval
pi = np.sqrt(e_pi)*z
zl = st.t.ppf(0.10, n-2) # z at 0.1
zu = st.t.ppf(0.90, n-2) # z at 0.9
ll = a + b*x1 + np.sqrt(e_pi)*zl # 10 %
ul = a + b*x1 + np.sqrt(e_pi)*zu # 90 %
In [42]:
import matplotlib.pyplot as plt
plt.plot(x,y,'ro', label='measured')
plt.plot(x1,ll,'--', label='10%')
plt.plot(x1,ul,'--', label='90%')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc='best')
plt.savefig('lin_regress.png')
In [ ]:
In [43]:
profile_time(1, 'k')
In [44]:
profile_time(4, 'k')
In [ ]:
In [ ]: