Add porosity to data projection

data/config.ini

@@ -2,6 +2,9 @@
 root=data
 base=Database_20220224.xyz
 hires=bathyhires.dat
+hstru=Hstru.dat
+poro=Poro.dat
+psize=Psize.dat
 hires_step=0.5
 
 [out]
@@ -11,6 +14,7 @@ out=bathy.npy
 step=1
 left=0
 right=150
+#plot=True
 
 [artha]
 lat=43.398450

data/processing/projection.py

@@ -5,6 +5,10 @@ import pathlib
 
 import numpy as np
 from scipy import interpolate
+try:
+    import matplotlib.pyplot as plt
+except ImportError:
+    plt = None
 
 from .lambert import Lambert
 
@@ -23,6 +27,9 @@ inp_root = pathlib.Path(config.get("inp", "root"))
 out_root = pathlib.Path(config.get("out", "root"))
 bathy_inp = out_root.joinpath(config.get("out", "sub"))
 hires_inp = inp_root.joinpath(config.get("inp", "hires"))
+hstru_inp = inp_root.joinpath(config.get("inp", "hstru"))
+poro_inp = inp_root.joinpath(config.get("inp", "poro"))
+psize_inp = inp_root.joinpath(config.get("inp", "psize"))
 bathy_out = inp_root.joinpath(config.get("out", "out"))
 
 log.info(f"Loading bathymetry from {bathy_inp}")
@@ -98,7 +105,37 @@ z[flt_x] = interpolate.griddata(
     (x[flt_x] - x[z_crossing] + bathy_hires[hires_crossing, 0]),
 )
 
+hstru_in = np.loadtxt(hstru_inp)[::-1]
+hstru = np.zeros(z.shape)
+hstru[flt_x] = interpolate.griddata(
+    (bathy_hires[:,0],),
+    hstru_in,
+    (x[flt_x] - x[z_crossing] + bathy_hires[hires_crossing, 0]),
+)
+
+poro_in = np.loadtxt(poro_inp)[::-1]
+poro = np.zeros(z.shape)
+poro[flt_x] = interpolate.griddata(
+    (bathy_hires[:,0],),
+    poro_in,
+    (x[flt_x] - x[z_crossing] + bathy_hires[hires_crossing, 0]),
+)
+
+psize_in = np.loadtxt(psize_inp)[::-1]
+psize = np.zeros(z.shape)
+psize[flt_x] = interpolate.griddata(
+    (bathy_hires[:,0],),
+    psize_in,
+    (x[flt_x] - x[z_crossing] + bathy_hires[hires_crossing, 0]),
+)
+
 np.savetxt(out_root.joinpath("bathy.dat"), z[::-1], newline=" ")
-np.savetxt(out_root.joinpath("hstru.dat"), np.zeros(z.shape), newline=" ")
+np.savetxt(out_root.joinpath("hstru.dat"), hstru[::-1], newline=" ")
-np.savetxt(out_root.joinpath("poro.dat"), np.zeros(z.shape), newline=" ")
+np.savetxt(out_root.joinpath("poro.dat"), poro[::-1], newline=" ")
-np.savetxt(out_root.joinpath("psize.dat"), np.zeros(z.shape), newline=" ")
+np.savetxt(out_root.joinpath("psize.dat"), psize[::-1], newline=" ")
+
+if plt is not None and config.getboolean("out", "plot", fallback=False):
+    fig, ax = plt.subplots()
+    ax.plot(-x, z, color="k")
+    ax.fill_between(-x, z+hstru, z, color="k", alpha=.2)
+    plt.show(block=True)