Trying to generate synthetic data with custom mesh

[jcmefra]

Problem description

Hello, I want the synthetic data generated with a Custom mesh but the example of ERT only shows how to use the custom mesh in the inversion, I want my synthetic model to be with a custom mesh as well. Thanks in advance.

Your environment

Please provide the output of print(pygimli.Report()) here. If that does not
work, please give provide some additional information on your:

Operating system: e.g. Windows, Linux or Mac?
Python version: e.g. 3.9, 3.10, etc.?
pyGIMLi version: Output of print(pygimli.__version__)
Way of installation: e.g. Conda package, manual compilation from source, etc.

--------------------------------------------------------------------------------
  Date: Thu Mar 14 22:48:33 2024 Hora est. Pacífico, Sudamérica

                OS : Windows
            CPU(s) : 8
           Machine : AMD64
      Architecture : 64bit
               RAM : 31.9 GiB
       Environment : Jupyter

  Python 3.9.18 | packaged by conda-forge | (main, Dec 23 2023, 16:29:04) [MSC
  v.1929 64 bit (AMD64)]

           pygimli : 1.4.6
            pgcore : 1.4.0
             numpy : 1.26.4
        matplotlib : 3.8.0
             scipy : 1.11.4
              tqdm : 4.66.2
           IPython : 8.18.1
           pyvista : 0.43.3

  Intel(R) oneAPI Math Kernel Library Version 2023.2-Product Build 20230613
  for Intel(R) 64 architecture applications
--------------------------------------------------------------------------------

Steps to reproduce

I'm building an iterable to create diverse random synthetic models:

import numpy as np
import pygimli as pg
import pygimli.meshtools as mt
from pygimli.physics import ert
import pandas as pd
import matplotlib.pyplot as plt
import random
import os

def generate_random_synth_data(index):
    base_path = 'simple_data'
    synth_images_path = os.path.join(base_path, 'synth_images')
    rhoa_distribution_path = os.path.join(base_path, 'rhoa_distribution')
    ertdata_path = os.path.join(base_path, 'raw')
    transformed_path = os.path.join(base_path, 'transformed')
    potential_distribution_path = os.path.join(base_path, 'potential_distribution')

    # Ensure directories exist
    os.makedirs(synth_images_path, exist_ok=True)
    os.makedirs(rhoa_distribution_path, exist_ok=True)
    os.makedirs(ertdata_path, exist_ok=True)
    os.makedirs(transformed_path, exist_ok=True)
    os.makedirs(potential_distribution_path, exist_ok=True)

    # Define the world randomly
    start_x = random.uniform(-50, -10)
    end_x = random.uniform(10, 50)
    depth = random.uniform(-30, -70)
    num_layers = random.randint(2, 5)

    world = mt.createWorld(start=[start_x, 0], end=[end_x, depth],
                       layers=[random.uniform(depth/3, 2*depth/3) for _ in range(num_layers - 1)], 
                       worldMarker=True)

    # Show the model and save the figure
    ax, _ = pg.show(world, showNodes=True)
    model_image_filename = f'ert_model_{index}.png'
    plt.savefig(os.path.join(synth_images_path, model_image_filename))
    plt.close()

    # Define electrode scheme randomly
    num_electrodes = 4
    elec_spacing = (end_x - start_x) / (num_electrodes - 1)
    elecs = np.linspace(start=start_x + elec_spacing, stop=end_x - elec_spacing, num=num_electrodes)
    schemes = ['slm']
    scheme = ert.createData(elecs=elecs, schemeName=random.choice(schemes))

    # Generate the mesh
    grid = pg.createGrid(x = np.linspace(start_x, end_x, 64), y = np.linspace(0, depth, 64))
    mesh = mt.createMesh(world, quality=34, grid=grid)

    # Define resistivity values randomly for layers
    default_resistivity = random.uniform(0.1, 1000)
    layer_resistivities = [random.uniform(0.1, 1000) for _ in range(num_layers)]
    rhomap = [[marker+1, res] for marker, res in enumerate(layer_resistivities)]
    rhomap.append([0, default_resistivity]) # Add default resistivity for any unspecified markers

    # Simulate measurements with random noise level
    noise_level = random.uniform(0, 0.05)  # Up to 5% noise
    noise_abs = random.uniform(0, 1e-4)  # Absolute noise
    data = ert.simulate(mesh, scheme=scheme, res=rhomap, noiseLevel=noise_level, noiseAbs=noise_abs)

    # Save simulation results
    data_filename = f'ertdata_{index}.dat'
    data.save(os.path.join(ertdata_path, data_filename))

    # Save the mesh with resistivity distribution as an image
    mesh_image_filename = f'ert_mesh_{index}.png'
    ax, _ = pg.show(mesh, data=rhomap, label=pg.unit('res'), showMesh=True)
    plt.savefig(os.path.join(synth_images_path, mesh_image_filename))
    plt.close()

    # Extracting node information for resistivity matrix
    node_coords = np.array([node.pos() for node in mesh.nodes()])
    x_coords = node_coords[:, 0]
    depths = node_coords[:, 1]
    resistivity_values = np.zeros(len(mesh.nodes()))
    
    for cell in mesh.cells():
        marker = cell.marker()
        resistivity = next(value for (mark, value) in rhomap if mark == marker)
        for node in cell.nodes():
            resistivity_values[node.id()] = resistivity

    # Combine into a matrix and save to CSV
    resistivity_matrix = pd.DataFrame({'X_coord': x_coords, 'Depth': depths, 'Resistivity': resistivity_values})
    resistivity_csv_filename = f'resistivity_distribution_{index}.csv'
    resistivity_matrix.to_csv(os.path.join(rhoa_distribution_path, resistivity_csv_filename), index=False)

    field = ert.simulate(mesh, scheme=scheme, res=rhomap, noiseLevel=noise_level, noiseAbs=noise_abs, sr=False, returnFields=True)
    pot = field[2]-field[1]

    ax, _ = pg.show(mesh, data=pot, cMap="RdBu_r",
                orientation='horizontal', label='Potential $u$',
                nCols=16, nLevs=9, showMesh=True, colorBar=False)

    cbar = plt.colorbar(ax.collections[0], ax=ax, orientation='vertical')

    #drawStreams(ax, mesh, pot, color='Black')

    plt.savefig(os.path.join(synth_images_path, f'ert_potential_{index}.png'))
    plt.close()

    # Extracting node information for potential matrix
    node_coords = np.array([node.pos() for node in mesh.nodes()])
    x_coords = node_coords[:, 0]
    depths = node_coords[:, 1]
    potential_values = pot

    # Combine into a matrix and save to CSV
    potential_matrix = pd.DataFrame({'X_coord': x_coords, 'Depth': depths, 'Potential': potential_values})
    potential_csv_filename = f'potential_distribution_{index}.csv'
    potential_matrix.to_csv(os.path.join(potential_distribution_path, potential_csv_filename), index=False)
    
# Generate datasets
n = 5  # Number of datasets to generate
for i in range(1, n+1):
    generate_random_synth_data(i)


for i in range(1, n+1):  # Adjust the range as needed
    # Define file paths
    base_path = 'simple_data/'
    dat_filename = f'{base_path}raw/ertdata_{i}.dat'
    mesh_image_filename = f'{base_path}synth_images/ert_mesh_{i}.png'
    model_image_filename = f'{base_path}synth_images/ert_model_{i}.png'
    resistivity_csv_filename = f'{base_path}rhoa_distribution/resistivity_distribution_{i}.csv'

    # Determine the number of electrodes from the file
    with open(dat_filename, 'r') as file:
        num_electrodes = int(file.readline().split()[0])

    # Load the electrode position data from the file, and compute the distance between electrodes
    electrode_position_data = pd.read_csv(dat_filename, sep='\t', skiprows=2, usecols=[0], nrows=num_electrodes, header=None)
    electrode_positions = electrode_position_data.values.flatten()
    distance_between_electrodes = electrode_positions[1] - electrode_positions[0]
    first_electrode_position = electrode_positions[0]

    # Determine the rows to skip
    rows_to_skip = 2 + num_electrodes + 2

    # Load the measurement data and convert electrode numbers to positions
    data = pd.read_csv(dat_filename, sep='\t', skiprows=rows_to_skip, names=['a', 'b', 'm', 'n', 'err', 'i', 'ip', 'iperr', 'k', 'r', 'rhoa', 'u', 'valid'])
    for electrode in ['a', 'b', 'm', 'n']:
        data[electrode] = (data[electrode] - 1) * distance_between_electrodes + first_electrode_position

    # Add a column for current (I) and fill it with 1.0 Amperes
    data['I'] = 1.0

    # Calculate the voltage (V = I * (rhoa / k)), which is basically the same as the resistance
    data['V'] = data['I'] * (data['rhoa'] / data['k'])

    # Select and rename columns for the new dataset, save the new dataset to a CSV file
    new_dataset = data[['a', 'b', 'm', 'n', 'I', 'V']]
    new_dataset.to_csv(f'{base_path}transformed/transformed_ertdata_{i}.csv', index=False)

    print(f"Processed dataset {i}")

print("Data generation and transformation complete.")
...

I think the part of interest is:

    # Generate the mesh
    grid = pg.createGrid(x = np.linspace(start_x, end_x, 64), y = np.linspace(0, depth, 64))
    mesh = mt.createMesh(world, quality=34, grid=grid)

Expected behavior

My synthetic models to have a 64x64 squared grid

Actual behavior

Grid is triangular:

If possible, please add one or more labels to your issue, e.g. if you expect that your issue is rather a question than a problem with the code, please add the label "question".

[halbmy]

I don't understand your mesh generation

grid = pg.createGrid(x=np.linspace(start_x, end_x, 64), 
                                  y=np.linspace(0, depth, 64))
mesh = mt.createMesh(world, quality=34, grid=grid)

The createMesh function does not have an argument grid, and it always creates a triangular mesh. You can directly use grid, or append a triangular boundary around it.

[jcmefra]

Thank you for the answer. I tried to just use grid but it outputs something like this:

I want to output something like the image I uploaded in the issue, but with a squared mesh so each node has a resistivity value and I can convert it to a matrix. Is that possible? Thank you!

[halbmy]

You mean you would like to populate the grid with markers (and then resistivities) like above?

for c in grid.cells():
    if c.center().y() > -20:
        c.setMarker(1)
    elif if c.center().y() > -30:
        c.setMarker(2)
...
pg.show(grid, markers=True)