As referenced in section 7.5.1 in Jackisch and Zehe (WRR, 2015), this is the testcase to reproduce a plot scale sprinkler experiment with the 2D version of echoRD - including a noised hydraulic conductivity definition.
Make sure to have numpy, pandas, scipy and echoRD installed.
In [1]:
import numpy as np
import pandas as pd
In [2]:
import scipy as sp
import matplotlib.pyplot as plt
import os, sys
In [3]:
#connect echoRD Tools
pathdir='../' #path to echoRD
lib_path = os.path.abspath(pathdir)
sys.path.append(lib_path)
import vG_conv as vG
from hydro_tools import plotparticles_t,hydroprofile,plotparticles_specht
In [4]:
obs=pd.read_csv('brspecht.dat',delimiter='\t')
obs.index=np.arange(-0.05,-1.,-0.1)
plt.figure(figsize=(18,6))
plt.subplot(131)
plt.imshow(obs.iloc[:,:10],cmap='Blues')
plt.colorbar()
plt.title('Br- Tracer Recovery\nProfile 1')
plt.subplot(132)
plt.imshow(obs.iloc[:,10:],cmap='Blues')
plt.title('Br- Tracer Recovery\nProfile 2')
plt.colorbar()
plt.subplot(133)
plt.plot(obs.iloc[:,:10].mean(axis=1),obs.index,label='Profile 1',c='b')
plt.errorbar(obs.iloc[:,:10].mean(axis=1),obs.index,xerr=obs.iloc[:,:10].std(axis=1),c='b')
plt.plot(obs.iloc[:,10:].mean(axis=1),obs.index,label='Profile 2',c='r')
plt.errorbar(obs.iloc[:,10:].mean(axis=1),obs.index,xerr=obs.iloc[:,10:].std(axis=1),c='r')
plt.legend(loc=4)
plt.xlabel('c(Br-)')
plt.ylabel('depth [m]')
plt.title('Br- Tracer Recovery [mg]\nLateral Mean with Var over Depth')
plt.tight_layout()
In [5]:
#connect to echoRD
import run_echoRD as rE
#connect and load project
[dr,mc,mcp,pdyn,cinf,vG]=rE.loadconnect(pathdir='../',mcinif='mcini_specht4y')
mcp.mcpick_out(mc,'specht4y_x2.pickle')
mc.advectref='Shipitalo'
[mc,particles,npart]=dr.particle_setup(mc)
precTS=pd.read_csv(mc.precf, sep=',',skiprows=3)
In [6]:
mc.prects=False
#theta=mc.zgrid[:,1]*0.+0.273
#[mc,particles,npart]=rE.particle_setup_obs(theta,mc,vG,dr,pdyn)
[thS,npart]=pdyn.gridupdate_thS(particles.lat,particles.z,mc)
[A,B]=plotparticles_t(particles,thS/100.,mc,vG,store=True)
In [7]:
mc.prects=False
precTS.tstart=60
precTS.tend=60+2.3*3600
precTS.total=0.02543
precTS.intense=precTS.total/(precTS.tend-precTS.tstart)
#DEBUG ONLY:
precTS.intense*=10.
precTS
Out[7]:
In [8]:
#NOISE KS
mu, sigma = 4.9e-6, 4.8e-6 # mean and standard deviation of observed ks at spechtacker
sl = np.random.normal(mu, 10**sigma, len(mc.soilgrid.ravel()))
subplot(121)
plot(mu*10**(sl),'.')
axhline(mu,c='r')
title('Noised Ks')
yscale('log')
subplot(122)
plot(10**(sl),'.')
axhline(1,c='r')
yscale('log')
title('multiplicator')
tight_layout()
In [9]:
ks_noise_factor=10**sl
#ks_noise_factor=1.
In [10]:
t_end=6.*3600.
saveDT=20
mc.FC[-1]=60
In [11]:
#1: MDA
#2: MED
infiltmeth='MDA'
#3: RWdiff
#4: Ediss
#exfiltmeth='RWdiff'
exfiltmeth='Ediss'
#5: film_uconst
#6: dynamic u
film=True
#7: maccoat1
#8: maccoat10
#9: maccoat100
macscale=1. #scale the macropore coating
clogswitch=False
In [ ]:
runname='echoRD_2D_Weiherbach_sprinkler_noise_'
drained=pd.DataFrame(np.array([]))
leftover=0
output=60. #mind to set also in TXstore.index definition
dummy=np.floor(t_end/output)
#prepare store arrays
TSstore=np.zeros((int(dummy),np.shape(thS)[0],np.shape(thS)[1]))
NPstore=np.zeros((int(dummy),len(mc.zgrid[:,1])+1))
colnames=['part1','part2','part3','part4','part5']
TXstore=pd.DataFrame(np.zeros((int(dummy),len(colnames))),columns=colnames)
t=0.
#loop through plot cycles
for i in np.arange(dummy.astype(int)):
#plot and store states
plotparticles_specht(particles,mc,pdyn,vG,runname,t,i,saving=True,relative=False)
[TXstore.iloc[i],NPstore[i,:]]=plotparticles_t(particles,thS/100.,mc,vG,runname,t,i,saving=True,store=True)
#store thS
TSstore[i,:,:]=thS
#check against br- profile
if (i==60) | (i==120) | (i==180) | (i==240) | (i==300):
z1=np.append(particles.loc[((particles.age>0.) & (particles.lat<=0.1)),'z'].values,obs.index)
z2=np.append(particles.loc[((particles.age>0.) & (particles.lat>0.1) & (particles.lat<=0.2)),'z'].values,obs.index)
z3=np.append(particles.loc[((particles.age>0.) & (particles.lat>0.2) & (particles.lat<=0.3)),'z'].values,obs.index)
advect_dummy1=(precTS.conc.values*mc.particleA)*(np.bincount(np.ceil((z1*-10).astype(float)).astype(np.int))[1:]-1)/mc.particleA
advect_dummy2=(precTS.conc.values*mc.particleA)*(np.bincount(np.ceil((z2*-10).astype(float)).astype(np.int))[1:]-1)/mc.particleA
advect_dummy3=(precTS.conc.values*mc.particleA)*(np.bincount(np.ceil((z3*-10).astype(float)).astype(np.int))[1:]-1)/mc.particleA
sim=pd.DataFrame([advect_dummy1,advect_dummy2,advect_dummy3],columns=obs.index).T
hydroprofile(obs,sim,fsize=(3, 6),xbound=[0.,0.03,3],ybound=[-1.,0.,5],ptitle='Bromid Recovery',saving=''.join(['./results/',runname,str(i),'BR.pdf']))
[particles,npart,thS,leftover,drained,t]=rE.CAOSpy_rundx_noise(i*output,(i+1)*output,mc,pdyn,cinf,precTS,particles,leftover,drained,6.,splitfac=4,prec_2D=False,maccoat=macscale,saveDT=saveDT,clogswitch=clogswitch,infilt_method=infiltmeth,exfilt_method=exfiltmeth,film=film,dynamic_pedo=True,ksnoise=ks_noise_factor)
In [ ]:
In [16]:
#pickle TSstore array - use only if a new model run was created (will overwrite file)
import cPickle as pickle
f = open(''.join([runname,'TS.pick']),"wb")
pickle.dump(pickle.dumps([pickle.dumps(TSstore),pickle.dumps(TXstore),pickle.dumps(NPstore)]), f, protocol=2)
f.close()
In [24]:
#in case there are different runs stored:
#runname='../resultsechoRD_2D_Weiherbach_1458z_'
#runname='NechoRD_2D_Weiherbach_1458A_'
In [25]:
#unpickle TSstore array
import cPickle as pickle
f = open(''.join(['./results/',runname,'TS.pick']),"rb")
[TS_store,TX_store,NP_store] = pickle.loads(pickle.load(f))
f.close()
TSstore=pickle.loads(TS_store)
TXstore=pickle.loads(TX_store)
NPstore=pickle.loads(NP_store)
In [26]:
np.shape(TSstore)
Out[26]:
In [39]:
import scipy.ndimage as spn
plt.figure(figsize=(18,6))
allplt=191
for i in np.arange(8):
k=allplt+i
plt.subplot(k)
plt.imshow(spn.filters.median_filter(TSstore[1+i**3,:,:],size=mc.smooth),vmin=0., vmax=100., cmap='Blues',origin='upper')
plt.title('t='+str(i**3)+'min',fontsize=16)
plt.subplot(k+1)
plt.colorbar()
plt.tight_layout()
plt.savefig(''.join(['./results/',runname,'moistdev.pdf']))
In [49]:
sim=pd.DataFrame(NPstore[-1,:]*precTS.conc.values*mc.particlemass*1000.,np.linspace(0.,mc.soildepth,len(NPstore[0,:])))
hydroprofile(obs,sim,fsize=(3, 6),xbound=[0.,0.025,3],ybound=[-1.,0.,5],ptitle='Bromid Recovery\nSimulation after 6h\nObservation after 24h\n',saving=''.join(['./results/',runname,str(360),'BR.pdf']))
In [47]:
def hydroprofile(obs,obs2=None,fsize=(6, 6),xbound=[0.,1.],ybound=[0.,1.],ptitle='Plot',saving=False,ccontrast=False,colors=None,compress=True):
'''
This is a rather simple function to plot hydrological data (profile) of a pandas data frame.
It is based on some excellent examples by Randy Olson and may need heavy adaption to your data.
(CC BY-NC-SA 4.0) jackisch@kit.edu
obs: dataframe to plot, index is y axis
mlab: labels for columns
mlabcols: columns to plot
fsize: Provide figure size as tuple.
saving: Provide a file name if you want it saving.
XXbound: Give bounds of axis and separations.
ptitle: Plot title
'''
# These are the "Tableau 20" colors as RGB.
tableau20 = [(31, 119, 180), (174, 199, 232), (255, 127, 14), (255, 187, 120),
(44, 160, 44), (152, 223, 138), (214, 39, 40), (255, 152, 150),
(148, 103, 189), (197, 176, 213), (140, 86, 75), (196, 156, 148),
(227, 119, 194), (247, 182, 210), (127, 127, 127), (199, 199, 199),
(188, 189, 34), (219, 219, 141), (23, 190, 207), (158, 218, 229)]
# Scale the RGB values to the [0, 1] range, which is the format matplotlib accepts.
for i in range(len(tableau20)):
r, g, b = tableau20[i]
tableau20[i] = (r / 255., g / 255., b / 255.)
if colors is None:
colors=tableau20
plt.figure(figsize=fsize)
# Remove the plot frame lines. They are unnecessary chartjunk.
ax = plt.subplot(111)
ax.spines["top"].set_visible(True)
ax.spines["bottom"].set_visible(False)
ax.spines["right"].set_visible(False)
ax.spines["left"].set_visible(True)
# Ensure that the axis ticks only show up on the bottom and left of the plot.
# Ticks on the right and top of the plot are generally unnecessary chartjunk.
ax.get_xaxis().tick_bottom()
ax.get_yaxis().tick_left()
# Limit the range of the plot to only where the data is.
# Avoid unnecessary whitespace.
plt.xlim(xbound[0], xbound[1])
plt.ylim(ybound[0], ybound[1])
# Make sure your axis ticks are large enough to be easily read.
# You don't want your viewers squinting to read your plot.
plt.yticks(np.linspace(ybound[0], ybound[1], ybound[2]), [str(x) for x in np.linspace(ybound[0], ybound[1], ybound[2])], fontsize=14)
plt.xticks(np.linspace(xbound[0], xbound[1], xbound[2]), [str(x) for x in np.linspace(xbound[0], xbound[1], xbound[2])],fontsize=14)
plt.ylabel('depth [m]')
plt.xlabel('c(Br-) [mg/l]')
ax.xaxis.grid(True)
ax.yaxis.grid(True)
# Now that the plot is prepared, it's time to actually plot the data!
# Note that I plotted the majors in order of the highest % in the final year.
if compress:
plt.plot(obs.mean(axis=1).values,obs.index, lw=2., color=colors[0], label='observed')
plt.errorbar(obs.mean(axis=1).values,obs.index,xerr=obs.std(axis=1).values, lw=1., color=colors[0])
offset=2
if (obs2 is not None):
plt.plot(obs2.mean(axis=1).values,obs2.index, lw=2., color=colors[2], label='simulated')
plt.errorbar(obs2.mean(axis=1).values,obs2.index,xerr=obs2.std(axis=1).values, lw=1., color=colors[2])
else:
for rank in range(len(obs.columns)):
# Plot each line separately with its own color, using the Tableau 20
# color set in order.
if ccontrast:
plt.plot(obs.iloc[:,rank].values,obs.index, lw=2., color=colors[rank*2], label=str(obs.columns[rank]))
else:
plt.plot(obs.iloc[:,rank].values,obs.index, lw=2., color=colors[rank], label=str(obs.columns[rank]))
offset=len(obs.columns)
if obs2 is not None:
for rank in range(len(obs2.columns)):
# Plot each line separately with its own color, using the Tableau 20
# color set in order.
if ccontrast:
plt.plot(obs2.iloc[:,rank].values,obs2.index, lw=2., color=colors[rank*2+offset], label=str(obs2.columns[rank]))
else:
plt.plot(obs2.iloc[:,rank].values,obs2.index, lw=2., color=colors[rank+offset], label=str(obs2.columns[rank]))
plt.title(ptitle,fontsize=16)
plt.legend(loc=4,frameon=False)
if saving!=False:
plt.savefig(saving, bbox_inches="tight")
In [ ]: