Simulate called with wrong resistivity array.

[nbillybgkv]

scheme = ert.createData(elecs=np.linspace(start=200, stop=300, num=21),
                           schemeName='dd')
data = ert.simulate(mesh, scheme=scheme, res=resistivity_values, noiseLevel=1,
                    noiseAbs=1e-6, seed=1337)

pg.info(np.linalg.norm(data['err']), np.linalg.norm(data['rhoa']))
pg.info('Simulated data', data)
pg.info('The data contains:', data.dataMap().keys())

pg.info('Simulated rhoa (min/max)', min(data['rhoa']), max(data['rhoa']))
pg.info('Selected data noise %(min/max)', min(data['err'])*100, max(data['err'])*100)

ERROR:BaseException Traceback (most recent call last)
Cell In[53], line 3
1 scheme = ert.createData(elecs=np.linspace(start=200, stop=300, num=21),
2 schemeName='dd')
----> 3 data = ert.simulate(mesh, scheme=scheme, res=resistivity_values, noiseLevel=1,
4 noiseAbs=1e-6, seed=1337)
6 pg.info(np.linalg.norm(data['err']), np.linalg.norm(data['rhoa']))
7 pg.info('Simulated data', data)

File ~.conda\envs\pg\Lib\site-packages\pygimli\physics\ert\ert.py:215, in simulate(mesh, scheme, res, **kwargs)
213 print(mesh)
214 print("res: ", res)
--> 215 raise BaseException(
216 "Simulate called with wrong resistivity array.")
218 if not isArrayData:
219 ret['rhoa'] = rhoa

BaseException: Simulate called with wrong resistivity array.

could you please tell me about this issue

[halbmy]

When opening a new issue, there is a template to be filled, e.g. specifying the version number etc. so that anyone can got a basis to help.

In your code, obviously resistivity_values is wrong. The output of the two print functions before the error should help.

[nbillybgkv]

Yes you are right..next time i will be careful/
The output is:
Mesh: Nodes: 238992 Cells: 237986 Boundaries: 476977
res: [130. 130. 130. ... 130. 130. 130.]

[halbmy]

So resistivity_values is a list? How long is it?

[nbillybgkv]

The Total number of resistivity values: 238992

[nbillybgkv]

it is numpy array.

[halbmy]

The length should match the cell number mesh.cellCount() and not the node number.

[nbillybgkv]

import numpy as np
1: first, I load a apparant resistivity image. Then, you define a range of resistivity values min_resistivity to max_resistivity. Next, calculate the minimum and maximum pixel values in the apparant resistivity image. and normalize them to the range [0, 1] using min-max normalization. Using these normalized pixel values, linearly scale them to the specified resistivity range, creating a resistivity image called resistivity_image. This process assigns resistivity values to each pixel in the segmented image based on its intensity,
1: # Define the range of resistivity values

min_resistivity = 42  # Ohm-meters
max_resistivity = 104  # Ohm-meters

#Calculate the range of pixel values in the apparent resistivity image
min_pixel_value = np.min(apparent resistivity_image)
max_pixel_value = np.max(apparent resistivity_image)

#Normalize the apparent resistivity image pixel values to the range [0, 1]
normalized_apparent resistivity_image = (apparent resistivty image - min_pixel_value) / (max_pixel_value - min_pixel_value)

#Map the Normalized pixel value to resistivity values within the specified range
resistivity_image = normalized_resistivity_image * (max_resistivity - min_resistivity) + min_resistivity

#Display the assigned resistivity values
print("Assigned resistivity values:")
print(resistivity_image).

output:Total number of resistivity values: 65536
Total number of pixels in the segmented image: 65536

Then in next step I Remove extra dimensions from the resistivity image

resistivity_image_2d = resistivity_image.squeeze()

#Check the shape of the 2D resistivity image
print("Shape of the 2D resistivity image:", resistivity_image_2d.shape)
total_resistivity_values = np.prod(resistivity_image_2d.shape)
print("Total number of resistivity values in the squeezed image:", total_resistivity_values)
output: Shape of the 2D resistivity image: (256, 256)
Total number of resistivity values in the squeezed image: 65536

 Then mesh is created based on the dimensions of a previously processed resistivity image. 
3: import numpy as np
import pygimli as pg
import pygimli.meshtools as mt

#Define mesh parameters
mesh_size = 1.0  # Size of mesh elements (in meters)

#Get dimensions of the resistivity image
x_dim, y_dim = resistivity_image_2d.shape

#Create a mesh grid based on dimensions and mesh size
x = np.arange(0, x_dim) * mesh_size
y = np.arange(0, y_dim) * mesh_size

#Create a mesh based on the mesh grid
mesh = mt.createMesh2D(x, y, markerType=1)

#Display the mesh
pg.show(mesh, label='Resistivity Mesh', showNodes=True, showMesh=True)
cell_count = mesh.cellCount()
#Print the cell count
print("Cell count:", cell_count)
output:Cell count: 65025

what is the issue that cell count 65025 is different than total number of resistivity values in the squeezed image: 65536.

[halbmy]

I have no idea what you're doing. Apparent resistivity is not a resistivity distribution but just a colored data table.

At any rate, for creating a mesh, you specify the edge positions. If you want to have 256 cells in each row, you need

x = np.arange(0, x_dim+1)

to make x and y 257 long and the mesh 256x256=65536. Otherwise you end up in 255x255=65025.

[nbillybgkv]

Actully i am trying to perform inversion through machine learning. For this purpose i am following the methdology of
Liu, B., Pang, Y., Jiang, P., Liu, Z., Liu, B., Zhang, Y., ... & Liu, J. (2023). Physics-driven deep learning inversion for direct current resistivity survey data. IEEE Transactions on Geoscience and Remote Sensing, 61, 1-11.