In de eerste plaats wordt een aantal libraries aangeroepen. Deze moeten zo nodig middels pip.install (library) worden geïnstalleerd. TIM en haar modules kunnen worden gevonden op Github. Het gehele script kan worden gevonden op deze Githubpagina
from concurrent.futures import ProcessPoolExecutor
import time
import os
import sys
if not sys.warnoptions:
import warnings
warnings.simplefilter("ignore")
import numpy as np
import pandas as pd
from ttim import *
from pylab import *
import shelve
import matplotlib_inline;
matplotlib_inline.backend_inline.set_matplotlib_formats('png', 'jpeg');
from datetime import datetime
Onderstaand blokje geeft de variabelen van de te tonen dataframes weer.
pd.set_option('display.max_columns', None)
pd.options.display.width=None
pd.set_option('display.max_rows', 10)
pd.set_option('display.expand_frame_repr', True)
pd.options.display.float_format = '{:,.2f}'.format
Lezen dict van k-waarden
Met behulp van de library shelve wordt de combinatie sedimentcode ingelezen als key en de doorlatendheid als value teruggegeven. Op deze wijze wordt de .gef boorbeschrijving omgevormd in een dataframe met k-waarden in [m/dag]. Een overzicht van de gehanteerde dictionary. Deze waarden zijn onder andere gebaseerd op ervaringswaarden, korrelgrootteanalyses, Pomper etc.
shelf = shelve.open("./Class", flag='r')
for key in shelf.keys():
Class = shelf[key]
shelf2 = shelve.open("./Class2", flag='r')
for key in shelf2.keys():
Class2 = shelf2[key]
shelf.close()
Lezen Gef-file
Hieronder staat het script om een .gef file uit het DINOloket te lezen en om te zetten in een dataframe met lagen van 1 meter dik. Bij dikkere lagen neemt de rekentijd af maar neemt de berekende doorlatendheid toe.
class GEF:
def __init__(self):
self._data_seperator = ' '
self._columns = {}
self.x = 0.
self.y = 0.
self.z = 0.
self.toplaag = []
self.ozlaag = []
self.typ = []
self.typ2 = []
self.ZFc = []
self.ZFres = []
def readFile(self, filename):
lines = open(filename, 'r').readlines()
for line in lines:
reading_header = True
for line in lines:
if reading_header:
self._parseHeaderLine(line)
else:
self._parseDataLine(line)
if line.find('#EOH') > -1:
if self._check_header():
reading_header = False
else:
return`
def _check_header(self):
if not 1 in self._columns:
return False
if not 2 in self._columns:
return False
return True
def _parseHeaderLine(self, line):
for xe in ['#COMMENT', 'Peil=', 'uitvoerder', 'materieel','WATERSTAND',
'opmerkingen#MEASUREMENTTEXT','==','= NAP','antropogeen']:
if xe in line:
return
if line.startswith(','):
return
if len(line.split()) == 0:
return
keyword, argline = line.split('=')
keyword = keyword.strip()
argline = argline.strip()
args = argline.split(',')
if keyword=='#XYID':
if float(args[1]) < 1e5:
args[1] = args[1]
else:
args[1]=args[1].replace('.','')
self.x = float(args[1])
args[2]=args[2].replace('.','')
self.y = float(args[2])
if (len(str(int(self.x))))>5:
self.x=int(self.x/pow(10,len(str(int(self.x)))-6))
if (len(str(int(self.y))))>5:
self.y=int(self.y/pow(10,len(str(int(self.y)))-6))
if self.x > 3e5:
self.x=self.x/10
elif keyword=='#ZID':
self.z = float(args[1])
elif keyword=='#COLUMNINFO':
column = int(args[0])
dtype = int(args[-1])
# if dtype==11:
# dtype = 10
self._columns[dtype] = column - 1
def _parseDataLine(self, line):
line = line.replace(' ','')
line = line.strip()
line = line.replace('|',' ')
line = line.replace(';',' ')
line = line.replace('!',' ')
line = line.replace(':',' ')
args = line.split()
for n, i in enumerate(args):
if i in ['9.9990e+003','-9999.9900','-1.0000e+007','-99999','-99999.0',
'-99999.00','-99999.000','-9.9990e+003','999.000', '-9999.99', '999.999',
'99','9.999','99.999', '999.9']:
args[n] = '0.1'
if len(line.split()) == 0:
return
zt = round(self.z-abs(float(args[self._columns[1]])),2)
zo = round(self.z-abs(float(args[self._columns[2]])),2)
tiep = args[9]
tiep = tiep.replace('\'','')
ZFres = np.nan
try:
ZF = args[11]
ZF = ZF.replace('\'','')
except IndexError:
ZF = np.nan
try:
ZFres = args[11]
ZFres = ZFres.replace('\'','')
except IndexError:
ZFres = np.nan
self.toplaag.append(zt)
self.ozlaag.append(zo)
self.typ.append(tiep)
self.ZFc.append(ZF)
Uiteindelijk zijn er 4 kolommen gemaakt met lagen van varierende dikte (conform boorbeschrijving) en met het bodemtype van die laag. Deze lagen zijn:
Hieronder worden de kolommen omgevormd tot een dataframe voor verder gebruik. In de def cijf staan twee correctie factoren f_PP en f_c , dit zijn handmatig aan te passen waarden (nu op 1 gesteld) om de doorlatendheid generiek te verhogen of te verlagen. Deze waarden worden gebruikt om met de uitkomst van de berekening de resulterende lijn op de punten te fitten. Deze correctie komt voort uit de persoon die de boorbeschrijving heeft gemaakt. Het bodemtype wat de een matig siltig zand noemt kan bij de nader licht kleiig fijn zand zijn, deze hebben een andere k-waarde in werkelijkheid. Deze correctiewaarde heeft mogelijk een afzettingsmilieu component. Zie het hoofdartikel.
def asNumpy(self):
return np.transpose(np.array([self.toplaag, self.ozlaag, self.typ, self.ZFc]))
def asDataFrame(self):
a = self.asNumpy()
return pd.DataFrame(data=a, columns=['tz','oz','laagsoort','ZFa'])
def cijf(self, filename):
df = self.asDataFrame()
df = df.sort_values('tz', ascending=False)
if df.empty:
return df
f_PP = 1
f_c = ff = 1
df['k_b'] = df['laagsoort'].map(Class)
df['zandF'] = df['ZFa'].map(Class2)
df['k'] = np.where(df['zandF']>0, df['k_b']*df['zandF']*f_PP, df['k_b']*f_PP)
df['tz'] = pd.to_numeric(df['tz'])
df['oz'] = pd.to_numeric(df['oz'])
df['c'] = 1/df['k']/f_c
df = df.sort_values('tz', ascending=False)
ymax=int((round(self.z,-1))+10)
ymin=ymax-150
df=df[df.tz >=ymin]
if df.empty:
return df
try:
dist1=round(df.iloc[3,0]-df.iloc[2,0],3)
dist2=round(df.iloc[-3,0]-df.iloc[-2,0],3)
dist = max(dist1, dist2)
dataq=np.array(np.arange(ymax,ymin,dist*-1))
except (IndexError, ValueError):
print(filename,'bestaat al')
return
dzend=int(round(df.loc[df.index[-1],'tz'],2))
Het dataframe wordt hieronder omgevormd tot een dataframe met lagen van 1 meter dik en de daarbij behorende gemiddelde k-waarde op basis van de boorbeschrijving. Het resulterende dataframe wordt geprint ter controle met de kolomnamen:
dataq=np.array(np.arange(ymax,dzend,-0.1))
dfq=pd.DataFrame(data=dataq, columns=['tz'])
dfnew=pd.merge(df,dfq,on='tz', how='outer')
dfnew=dfnew.sort_values(by=['tz'], ascending=False)
dfnew.iloc[[0], [1]] = -1
dfnew['k'] = dfnew['k'].ffill()
dfnew['k'] = pd.to_numeric(dfnew['k'])
df = dfnew
df=df.reset_index(drop=True)
Tfa = []
Val = []
Z = int((round(self.z,-1))+10)
df2 = df.set_index("tz", drop =True)
for Tf in range(Z, dzend, -1):
df3 = df2.loc[Tf:(Tf-1):, 'k'].mean()
Tfa.append(Tf)
Val.append(df3)
a=pd.DataFrame(Tfa, columns=['Top_zand'])
a['OZ_zand']=a['Top_zand']-1
a['k'] = Val
a=a.dropna()
print(a)
Data
In dit deel van het script worden de parameters van de pompproef hard ingevoerd, daarnaast is er nog de parameter Cf, dit is een vermenigvuldigingswaarde om de curve aan de hand van de berging te fitten, en de Anif, hetzelfde maar dan voor de anisotropie. Deze twee factoren zijn vaak de hier aangegeven waarden. Bij een unconfined WVP kan deze waard soms op 50 liggen waarmee de berging het porien volume benadert. De filterquote is de lengte van het onttrekkingsfilter in relatie tot de dikte van het WVP (de onvolkomendheid van het filter)
t = (3) # Lengte pompproef in dagen
print('duur pompproef ' , round(t,1), ' dagen')
Qm1 = 24*39 # Debiet dag
Tf1 = -29 # Top filter in m NAP zie gebiedsdossier
Bf1 = -37 # Onderzijde filter in m NAP
fT = self.z # Onderzijde freatisch niveau (GLG?))
Cf = 0.2 # kleiner is lagere S
Anif = 0.5 # hoger is hogere anisotropie
TWVP = -28 # Top Watervoerend pakket in [m NAP]
OWVP = -39 # Onderzijde Watervoerend pakket in [m NAP]
print(filename)
print('Filterquote = ', round(((abs(Bf1-Tf1))/(abs(OWVP-TWVP))),2))
De code hieronder creëert een aantal waarden per modellaag. Aan de hand van een tweetal algoritmes (zie artikel) voor de Berging (S) per meter (=Ss) en de Anisotropie per meter (Kzkh) worden aan de dataframe kolommen toegevoegd. Tenslotte wordt het resulterende dataframe geprint. Dit wordt gebruikt als inputmodel voor de TTIM module
a['Anif'] = Anif
a=a.drop(columns=['OZ_zand'])
a.columns=['Top', 'k', 'Anif']
a=a.dropna()
a['k'] = np.where(a['k']<1, a['k']/ff, a['k'])
a['kzkh'] = 2/a['Anif']/(33.653*np.exp(-0.066*(a['k'])))
a['kzkh'] = np.where(a['kzkh']> 1, 1, a[ 'kzkh'])
a['an'] = (np.round(1/a['kzkh'],2)).astype('float')
a['S'] = (a['k']/(24*60*60))*a['kzkh']*Cf
a['S'] = np.where(a['Top'] >fT, 0.05, a['S'])
a['S'] = a['S'].apply(lambda x: '%.1e' % x).values.tolist()
S = a['S'].tolist()
S.append(1e-5)
a = a.reset_index(drop=True)
a['kd'] = a['d']*a['k']
a['c'] = 1/a['k']
print(a)
Om TTIM correct te laten werken worden numpy arrays gemaakt. Hierbij is de hoogteligging 1 element langer dan bijvoorbeeld de k-waarde.
Z= a['Top'].tolist()
Z.append(dzend-1)
Z.append(dzend-2)
Kaq= a['k'].tolist()
Kaq.append(1e-3)
kzkh=a['kzkh'].tolist()
kzkh.append(0.15)
tmin=1e-4
tmax=1e10
De code hieronder bepaalt de modellagen waar de filter in is geplaatst. Verder wordt het deel waar het WVP inzit als los dataframe geknipt voor de bepaling van de KD waarde van het WVP en de totale Berging (S) van het WVP. Alle dataframe waarden wordne naar een aparte xls file gekopieerd en in de werkdirectory opgeslagen.
b=a
d=a
b=b.drop(b[(b.Top > Tf1)].index)
b=b.drop(b[(Bf1 >= b.Top)].index)
d =d.drop(d[(d.Top > TWVP)].index)
d =d.drop(d[(OWVP >= d.Top)].index)
Ld=list(d.index.values) #WVP_lagen
d['S'] = pd.to_numeric(d['S'])
St = d['S'].sum()
if St>0.35:
St=0.35
c1t = df.drop(df[(df.tz > TWVP)].index)
c1b = c1t.drop(df[(df.tz < OWVP)].index)
KDl1=list(c1b.index.values)
dikte1 = abs(np.round(float(c1b.tz.min()-c1b.tz.max()),1))
print('Dikte WVP [m] = ', dikte1)
Berging = "{:.1E}".format(St)
St_m = d['S'].sum()/dikte1
Berging_m = "{:.1E}".format(St_m)
print('Berging [-] = ', Berging)
print('Berging [1/m] = ', Berging_m)
print('Uurdebiet ', int(np.round(Qm1/24,1)), '[m3/uur] Dagdebiet', int(np.round(Qm1,1)),' [m3/dag]')
ab = (a[a.Top < TWVP])
ac = (ab[ab.Top > OWVP])
KD = int(ac.kd.sum())
KD1= round(KD/dikte1,1)
print('KD WVP [m2/dag] = ', KD)
print('k_gem [m/dag] = ', KD1)
Lb=list(b.index.values) #Pomplagen1
print('Filterlagen1 = ', Lb)
a['c'] = a['c'].astype('int')
a.to_excel('invoer.xlsx')
az= a[a.k>1.0]
az = az[az['Top'].between(OWVP, TWVP)]
print('Berekende anisotropie zand: ', round(az.an.median(),1))
print('k_factor = ', f_PP)
print('fijnfactor = ', ff)
In het codeblok hieronder wordt de TTIM module aangeroepen en gevuld met de rekenwaarden. Een library wordt gebruikt om de berekeningen door middel van parallel rekenen te versnellen. Uiteindelijk wordt in de regel met ml.solve() de daadwerkelijke som gemaakt.
ml=Model3D(kaq=Kaq, z=Z, Saq= S, kzoverkh=kzkh, phreatictop=False,
tmin=tmin, tmax=tmax, M=20)
w1 = Well(ml, xw=0, yw=0, rw=0.30, tsandQ=[(0, Qm1)], layers=Lb)
aantal_cores = os.cpu_count()
pool = ProcessPoolExecutor(max_workers=aantal_cores)
list(pool.map(ml.solve()))
Voor de website wordt een .md file gemaakt van de laagopbouw en parameters.
a['Top'] = np.where(a['Top']==int(0), 'NAP', a['Top'])
mdtable = a[['Top','k', 'an', 'S', 'kd', 'c']].to_markdown(stralign='right',
floatfmt=(".0f",".0f", ".0f",".1f",".1e", ".0f",".0f",".0f" ))
mdtable = mdtable.replace(' 0', ' ')
with open('invoer.md', 'w') as file:
file.write(mdtable)
Tenslotte wordt de uitkomst van de berekening gevisualiseerd en wordt de fit met de meetwaarden vergeleken.
laag = [Lb[1]]
za=ml.plots.head_along_line(x1=0.1, x2=2500, y1=0, y2=0, npoints=10000, t=t,
layers=[laag[0]], lw=2, grid=True)
Bronstand = round(w1.headinside(t).max(),2)
print ('Waterpeil in PP = ', Bronstand)
leg=plt.legend(round(a.loc[a.index[laag], 'Top'],1), loc=(0.85,0.05),
fontsize = 8, title='Filter op [m NAP]')
plt.semilogx()
plt.xlim(0.1,2500)
plt.ylim(-8,0)
plt.ylabel('Verlaging [m]')
plt.xlabel('Afstand tot bron [m]')
title=('KD berekend= '+str(int(KD))+' [m2/dag]')
plt.title(title, fontsize=12, pad=10)
plt.plot(0.4, -6.40, 'b+', markersize=10)
plt.plot(22, -2.43, 'b+', markersize=10)
plt.plot(40, -1.90, 'b+', markersize=10)
plt.plot(80, -1.38, 'b+', markersize=10)
plt.plot(150, -0.88, 'b+', markersize=10)
plt.plot(255, -0.57, 'b+', markersize=10)
plt.plot(0.12, Bronstand,'ro', markersize=5)
plt.savefig('test.png', bbox_inches='tight')
plt.show()
plt.close()
az= a[a.k>1.0]
az = az[az['Top'].between(OWVP, TWVP)]
print('Berekende anisotropie zand: ', round(az.an.median(),1))
end = time.time()
hours, rem = divmod(end-strt, 3600)
minutes, seconds = divmod(rem, 60)
print('Rekentijd: '+"{:0>2}:{:0>2}:{:0>2}".format(int(hours),int(minutes),int(seconds)))
print('Start berekening: ', datetime.now().strftime('%Y-%m-%d %H:%M:%S'))
strt=time.time()
for filename in os.listdir(os.getcwd()):
if filename.endswith ('.GEF') or filename.endswith ('.gef'):
if __name__=="__main__":
g=GEF()
g.readFile(filename)
g.cijf(filename)
plt.show()
plt.close()