Map porepressure from a csv file with only positive values to FLAC3D gridpoints.

map_pp_positive.py

import itasca as it
from itasca import gpa
import numpy as np
from scipy.constants import foot
from collections import defaultdict
import os
from scipy.spatial import cKDTree
import scipy

if os.path.exists("data.npy"):
    data = np.load("data.npy")
else:
    data = np.loadtxt("My_pore_pressure_data.csv",
                      delimiter=",", usecols=(1,2,3,4,5), skiprows=1)
    np.save("data.npy", data)

# write for Paraview so we can visualize
data *= [foot,foot,foot,1,1000*9.81*foot] # convert to m,m,m,s,Pa
X,Y,Z,time,head = data.T



f3_gridpoints = gpa.pos()


def draw_fragment(gname, points, faces):
    f = open(gname, "w")
    print("ITASCA GEOMETRY3D", file=f)
    print("NODES", file=f)
    for j, p in enumerate(points):
        x,y,z = p
        print(j+1, x, y, z, "EXTRA 1 0", file=f)

    print("POLYS", file=f)
    for j, t in enumerate(faces):
        n0, n1, n2 = t
        print(j+1, "NODE", n0+1, n1+1, n2+1, "EXTRA 1 0", file=f)
    f.close()


# this type of data has repeating x,y coordinates
# we bin them first


times = list(set(time))
for t in times:
    mask = time==t
    print("running t=", t)
    gridded = defaultdict(list)
    for row in data[mask]:
        x,y,z,_,head = row
        gridded[(x,y)].append((z,head))

    new_points, new_head = [], []
    free_surface = []
    for x_y, z_head in gridded.items():
        z_head = list(set(z_head))
        z_head.sort(key=lambda _:_[0])
        z0, h0 = z_head[-2]
        z1, h1 = z_head[-1]
        assert not z1==z0
        assert not h1==h0
        # look at the vertical pore pressure gradient and project upwards to find the pressure = zero height
        slope = (h1-h0)/(z1-z0)
        # we can calculate the slope or use the (known) hydrostatic gradient
        slope = -1_000*9.81
        assert slope < 0
        delta_h = -h1/slope
        new_points.append((*x_y, z1 + delta_h*2))
        # because we are using nearest neighbor, we need to double delta h so we define the zero surface correctly.
        free_surface.append((*x_y, z1 + delta_h))
        new_head.append(0)
    all_points = np.vstack((data[mask][:,:3], np.array(new_points)))
    all_head = np.hstack((data[mask][:,-1], np.array(new_head)))
    #all_head
    distances, indices = cKDTree(all_points).query(f3_gridpoints)
    print(np.histogram(distances))
    gp_heads = all_head[indices]
    print(f"saving pp_t_{int(t)}.npy")
    np.save(f"gp_pp_t_{int(t)}.npy", gp_heads)
    new_points = np.array(new_points)
    mesh = scipy.spatial.Delaunay(np.array(free_surface)[:,:2])
    draw_fragment(f"free_surface_t_{int(t)}.geom", np.array(free_surface), mesh.simplices)

This is for a special case of porepressure data that is kind of common. All the points are aligned on a repeating x,y grid but they can be at any height. There are no zero values so the nearest neighbor interpolation does not directly work. The procedure is to identify all the columns of points with the same x and y coordinates, sort them by height, look at the highest one, extrapolate upwards to find the phreatic surface, and add those as zero pressure points to the data. Then we do the nearest neighbor interpolation.