diff --git a/README.md b/README.md
index 9ced9806b4ca8d216f5d067c97a993a374a27476..8949acf72b1a737a9733a28b0a882c2e18fce7b4 100644
--- a/README.md
+++ b/README.md
@@ -1,3 +1,67 @@
# adists_sediment_cores
-Python script to extract numerical sediment cores from adists binary results
\ No newline at end of file
+Python script to extract numerical sediment cores from adists binary results
+
+## données:
+
+* h[0:tmax] : liste des hauteurs totales de sédiment à chaque pas de temps
+* masse[0:nb_sediment*nb_affluents][0:tmax] : la masse de chaqun des sédiments qui constituent le fond a chaque pas de temps
+adis-ts n'a qu'une seule couche et ne mélange pas les sédiments : on a besoin d'un algo pour recréer ces couches
+note : on part du principe qu'on n'a rien pour relier la masse à la hauteur
+
+## sorties :
+
+* la hauteur de chaque couche
+* sa proportion en chaque sédiment
+
+## algorythm : (pseudo-code, ne fonctionne pas)
+
+```python
+layers = [] # liste de couches, vide pour l'instant
+l = 0 # numéro de la couche courante
+for t in range(nb_timesteps):
+ # on crée une nouvelle couche en fin de liste
+ # une couche contient:
+ # [0] : la hauteur de sédiments
+ # [1] : la masse des sédiments dans la couche
+ # [2:nb_sediments+2] : la proportion de chaque classe de sédiments dans la couche
+ layers.append(np.zeros(nb_sediments+2,float))
+ layers[l][0] = h[t]-h[t-1] # hauteur potentiellement négative
+ layers[l][1] = 0.0
+ for s in range(nb_sediments):
+ # pour l'instant on met la masse de chaque sédiment dans layers[l][2:nb_sediments+2]
+ # la masse d'une classe de sédiment est la somme des masses des classes correspondantes de chaque affluent
+ # il s'agit d'une diférence avec le pas de temps précédent, ces masses peuvent donc être négatives si il y a érosion
+ layers[l][s+2] = sum(masse[s+j][t]-masse[s+j][t-1] for j in range(0,nb_afluent*nb_sediments,nb_sediments)
+ residual = layers[l][s+2]
+ if residual < 0.0:
+ # il y a eu érosion pour cette classe
+ # dans notre carotte on ne peut pas avoir une masse négative
+ layers[l][s+2] = 0.0
+ # on ajoute la masse de cette classe à la masse totale de la couche
+ layers[l][1] += layers[l][s+2]
+ l2=l
+ # si residual est négatif alors on doit aller chercher du sédiment dans les couches du dessous
+ while residual < 0.0:
+ l2 -= 1
+ # calcul préliminaire
+ diff = min(-1*residual, layers[l2][s+2]) # diff > 0 : ce que l'on doit prélever
+ ratio = (layers[l2][1]-diff)/layers[l2][1] # ratio < 1 : pour réduire l'épaisseur de la couche.
+ dh = layers[l2][0]-(ratio*layers[l2][0]) # dh > 0
+ # modif de la couche
+ residual += layers[l2][s+2]
+ layers[l2][s+2] -= diff
+ layers[l2][1] -= diff
+ layers[l2][0] -= dh
+ layers[l][0] += dh # doit finir par compenser les hauteurs négatives
+ # fin du while sur residual
+ # fin de la boucle sur les classes de sédiments
+ l += 1
+# on a fini, on va refaire une boucle sur les couches pour obtenir les proportions
+for layer in layers:
+ for s in range(nb_sediments):
+ if layers[l][1] > 0.0: # normalement jamais négatif
+ # on transforme les masses en proportions
+ layer[s+2] /= layer[1]
+```
+
diff --git a/sed_cores/options.py b/sed_cores/options.py
index 92bbef055ed1481ffe99bf5c02d01cf3748e3781..b78b9618ae7fb673cef7a638f9f22548413f3ef6 100644
--- a/sed_cores/options.py
+++ b/sed_cores/options.py
@@ -2,6 +2,9 @@ extract_adists_data = True
adists_bin_data_dir = "."
adists_extracted_data_dir = "./adists_out"
adists_exe = "adists"
+plot_instant_accumulation = False
+plot_positive_accumulation = False
+plot_negative_accumulation = False
output_accumulation_to_csv = True
output_accumulation_to_png = True
output_cores_to_csv = True
diff --git a/sed_cores/sediment_cores.py b/sed_cores/sediment_cores.py
index 07e6bc73e97631273c7b0a56a09934a158f236b8..70753cc9df7c8fbae2095e5e5e32a661e32846ce 100644
--- a/sed_cores/sediment_cores.py
+++ b/sed_cores/sediment_cores.py
@@ -127,6 +127,9 @@ from options import extract_adists_data
from options import adists_bin_data_dir
from options import adists_extracted_data_dir
from options import adists_exe
+from options import plot_instant_accumulation
+from options import plot_positive_accumulation
+from options import plot_negative_accumulation
from options import output_accumulation_to_csv
from options import output_accumulation_to_png
from options import output_cores_to_csv
@@ -169,14 +172,12 @@ if extract_adists_data:
os.mkdir(adists_extracted_data_dir)
os.chdir(adists_extracted_data_dir)
print("Extracting data from Adis-TS to CSV files...")
- print("we are in",os.getcwd())
output = ""
if not verbose: output = " > /dev/null"
for bed in list_beds:
for file in list_files:
for i,pk in enumerate(list_pk):
- print("Extracting data from "+adists_bin_data_dir+"/"+file+".bin")
os.system( adists+" -b "+adists_bin_data_dir+"/"+file+".bin "+bed+"dt "+str(pk)+" "+file+"_"+bed+str(i)+output)
os.chdir(cwd)
@@ -232,33 +233,42 @@ if (not os.path.isdir("./accumulation")):
os.chdir("./accumulation")
values_accumulated = dict()
-values_accumulated_tot = dict()
+values_accumulated_pos = dict()
+values_accumulated_neg = dict()
nb_layers = dict()
for bed in list_beds:
if verbose: print(" Bed: "+bed)
if def_numpy:
nb_layers[bed] = np.zeros(nb_pk,int)
- values_accumulated[bed] = np.zeros((nb_pk,nb_timesteps),float)
- values_accumulated_tot[bed] = np.zeros((nb_pk,nb_timesteps),float)
+ values_accumulated[bed] = np.zeros((nb_files,nb_pk,nb_timesteps),float)
+ values_accumulated_pos[bed] = np.zeros((nb_pk,nb_timesteps),float)
+ values_accumulated_neg[bed] = np.zeros((nb_pk,nb_timesteps),float)
else:
nb_layers[bed] = create_int_list([nb_pk])
- values_accumulated[bed] = create_float_list([nb_pk,nb_timesteps])
- values_accumulated_tot[bed] = create_float_list([nb_pk,nb_timesteps])
+ values_accumulated[bed] = create_float_list([nb_files,nb_pk,nb_timesteps])
+ values_accumulated_pos[bed] = create_float_list([nb_pk,nb_timesteps])
+ values_accumulated_neg[bed] = create_float_list([nb_pk,nb_timesteps])
val_acc = values_accumulated.get(bed)
- val_acc_tot = values_accumulated_tot.get(bed)
+ val_acc_pos = values_accumulated_pos.get(bed)
+ val_acc_neg = values_accumulated_neg.get(bed)
val = values.get(bed)
layers = nb_layers.get(bed)
for pk in range(nb_pk):
if verbose: print(" PK: "+str(pk))
- val_acc[pk][0] = val[0][pk][0]
- val_acc_tot[pk][0] = val[0][pk][0]
+ for j in range(nb_files):
+ val_acc[j][pk][0] = val[j][pk][0]
+ val_acc_pos[pk][0] = val[0][pk][0]
+ val_acc_neg[pk][0] = val[0][pk][0]
for t in range(nb_timesteps-1):
- val_acc[pk][t+1] = val[0][pk][t+1] - val[0][pk][t]
- if val[0][pk][t+1]>val[0][pk][t]:
+ for j in range(nb_files):
+ val_acc[j][pk][t+1] = val[j][pk][t+1] - val[j][pk][t]
+ if val_acc[0][pk][t+1]>0.0:
layers[pk] += 1
- val_acc_tot[pk][t+1] = val_acc_tot[pk][t] + val_acc[pk][t+1]
+ val_acc_pos[pk][t+1] = val_acc_pos[pk][t] + val_acc[0][pk][t+1]
+ val_acc_neg[pk][t+1] = val_acc_neg[pk][t]
else:
- val_acc_tot[pk][t+1] = val_acc_tot[pk][t]
+ val_acc_pos[pk][t+1] = val_acc_pos[pk][t]
+ val_acc_neg[pk][t+1] = val_acc_neg[pk][t] + val_acc[0][pk][t+1]
os.chdir(cwd)
if verbose: print(" Computing accumulation done")
@@ -271,27 +281,29 @@ if output_accumulation_to_csv:
os.chdir("./accumulation")
if def_numpy:
- row_data=np.zeros((nb_timesteps,4),float)
+ row_data=np.zeros((nb_timesteps,5),float)
else:
- row_data = create_float_list([nb_timesteps,4])
+ row_data = create_float_list([nb_timesteps,5])
for i in range(nb_timesteps):
row_data[i][0]=list_timesteps[i]
for bed in list_beds:
val_acc = values_accumulated.get(bed)
- val_acc_tot = values_accumulated_tot.get(bed)
+ val_acc_pos = values_accumulated_pos.get(bed)
+ val_acc_neg = values_accumulated_neg.get(bed)
val = values.get(bed)
layers = nb_layers.get(bed)
for pk in range(nb_pk):
for t in range(nb_timesteps):
- row_data[t][1]=val_acc[pk][t]
- row_data[t][2]=val_acc_tot[pk][t]
- row_data[t][3]=val[0][pk][t]
+ row_data[t][1]=val_acc[0][pk][t]
+ row_data[t][2]=val_acc_pos[pk][t]
+ row_data[t][3]=val_acc_neg[pk][t]
+ row_data[t][4]=val[0][pk][t]
if layers[pk]>0:
if verbose: print(" Writing values_accumulated_"+list_pk[pk]+"_"+bed+".csv")
with open("values_accumulated_"+list_pk[pk]+"_"+bed+".csv", "w") as f:
writer = csv.writer(f)
- writer.writerow(["timestep","instant_accumulation","positive accumulation","total_accumulation"])
+ writer.writerow(["timestep","instant_accumulation","positive accumulation","negative accumulation","total_accumulation"])
for t in range(nb_timesteps):
writer.writerow(row_data[t][:])
f.close()
@@ -306,7 +318,8 @@ if def_matplotlib and output_accumulation_to_png:
os.chdir("./accumulation")
for bed in list_beds:
val_acc = values_accumulated.get(bed)
- val_acc_tot = values_accumulated_tot.get(bed)
+ val_acc_pos = values_accumulated_pos.get(bed)
+ val_acc_neg = values_accumulated_neg.get(bed)
layers = nb_layers.get(bed)
val = values.get(bed)
for pk in range(nb_pk):
@@ -317,9 +330,13 @@ if def_matplotlib and output_accumulation_to_png:
ax.set_xlabel("time")
ax.set_ylabel("height")
#ax.set_title("Sediment accumulated at pk"+list_pk[pk]+" in bed "+bed)
- plt.plot(list_timesteps[:],val_acc[pk][:],label="instant accumulation")
- plt.plot(list_timesteps[:],val_acc_tot[pk][:],label="positive accumulation")
- plt.plot(list_timesteps[:],val[0][pk][:],label="total accumulation")
+ if plot_instant_accumulation:
+ plt.plot(list_timesteps[:],val_acc[0][pk][:],label="instant accumulation")
+ if plot_positive_accumulation:
+ plt.plot(list_timesteps[:],val_acc_pos[pk][:],label="positive accumulation")
+ if plot_negative_accumulation:
+ plt.plot(list_timesteps[:],val_acc_neg[pk][:],label="negative accumulation")
+ plt.plot(list_timesteps[:],val[0][pk][:],label="accumulation")
plt.legend()
plt.title(label="Sediment accumulated at pk"+list_pk[pk]+" in bed "+bed+"\n", fontsize = 20)
@@ -343,33 +360,54 @@ if def_numpy:
val_temp1 = np.zeros(nb_sediments,float)
else:
val_temp1 = create_float_list([nb_sediments])
-first_row=["height","layer_size"]
+first_row=["height","mass"]
for sed in list_sediments : first_row.append(sed)
for bed in list_beds:
if verbose: print(" Bed: "+bed)
layers[bed] = []
val_acc = values_accumulated.get(bed)
val = values.get(bed)
- val_acc_tot = values_accumulated_tot.get(bed)
+ val_acc_pos = values_accumulated_pos.get(bed)
for pk in range(nb_pk):
if verbose: print(" PK: "+str(pk))
layers[bed].append([])
l = 0
for t in range(nb_timesteps):
- if val_acc[pk][t] > 0.0:
- if def_numpy:
- layers[bed][pk].append(np.zeros(nb_sediments+2,float))
- else: #def_numpy
- layers[bed][pk].append(create_float_list([nb_sediments+2]))
- layers[bed][pk][l][0] = val_acc[pk][t]
- layers[bed][pk][l][1] = val_acc_tot[pk][t]
- #print(layers[bed][pk][l])
+ if def_numpy:
+ layers[bed][pk].append(np.zeros(nb_sediments+2,float))
+ else: #def_numpy
+ layers[bed][pk].append(create_float_list([nb_sediments+2]))
+ layers[bed][pk][l][0] = val_acc[0][pk][t] # potentiellement négative
+ layers[bed][pk][l][1] = 0.0
+ for s in range(nb_sediments):
+ layers[bed][pk][l][s+2] = sum(val_acc[s+j][pk][t] for j in range(1,nb_files,nb_sediments))
+ residual = layers[bed][pk][l][s+2]
+ if residual < 0.0:
+ layers[bed][pk][l][s+2] = 0.0
+ layers[bed][pk][l][1] += layers[bed][pk][l][s+2]
+ l2=l
+ while residual < 0.0:
+ l2 -= 1
+ # calcul préliminaire
+ diff = min(-1*residual, layers[bed][pk][l2][s+2]) # diff > 0
+ ratio = (layers[bed][pk][l2][1]-diff)/layers[bed][pk][l2][1] # ratio < 1
+ dh = layers[bed][pk][l2][0]-(ratio*layers[bed][pk][l2][0]) # dh > 0
+ # modif de la couche
+ residual += layers[bed][pk][l2][s+2]
+ layers[bed][pk][l2][s+2] -= diff
+ layers[bed][pk][l2][1] -= diff
+ layers[bed][pk][l2][0] -= dh
+ layers[bed][pk][l][0] += dh
+ if layers[bed][pk][l][1] <= 0.0:
+ layers[bed][pk].pop() # pas besoin d'une couche vide
+ else:
l += 1
- val_temp2 = sum(val[j][pk][t] for j in range(1,nb_files))
+ for layer in layers[bed][pk]:
+ if layers[1] > 0.0:
for s in range(nb_sediments):
- val_temp1[s] = sum(val[s+j][pk][t] for j in range(1,nb_files,nb_sediments))
- if val_temp2>0.0 and val_temp1[s]>0.0 :
- layers[bed][pk][l-1][s+2] = val_temp1[s]/val_temp2
+ # on transforme les masses en proportions pour l'affichage
+ layer[s+2] /= layer[1]
+
os.chdir(cwd)
if verbose: print("Computing sediment cores done")