DC-2d inversion with Wenner

Hi, I’m a geophysics graduate student practicing doing DC 2D inversion, I’m trying to use the code example provided on your website 2.5D DC Resistivity Least-Squares Inversion — SimPEG 0.14.3 documentation to do some SIMPEG inversion.

However, when I try to produce the Pesudosection of the data by using dc.utils.plot_pseudoSection(), it always jumps the error message:
IndexError: index 0 is out of bounds for axis 1 with size 0
Or
sometimes gives this: ValueError: all the input array dimensions for the concatenation axis must match exactly, but along dimension 0, the array at index 0 has size 573 and the array at index 1 has size 1

I have checked as many times as I can, I couldn’t find the reason for this error.
I’m trying to use Wenner Geometry instead of dipole-dipole.
If you could provide an example of how to do Wenner for 2D inversion, it would be helpful.

Thanks in advance.

Hi Yun,
if you could share your code that would be easier to resolve.

Hi Thibaut, here is the code. the observation data format is like below(I can’t upload the data from here).

10m multinodes wenner
10.00
1
573
1
0
15 10 31.4
30 20 31.0
45 30 30.1
60 40 27.6
75 50 24.8
90 60 22.5
105 70 21.9
120 80 19.4
135 90 16.8
25 10 31.2
40 20 31.6
35 10 29.7
50 20 30.8
65 30 29.9
80 40 29.1
95 50 26.0
155 90 18.6
45 10 29.0
60 20 29.8
120 60 22.1
55 10 27.8
70 20 29.0
85 30 29.9
100 40 28.9
.
.
.
If you need this data from me, let me know, I’ll send it to you through other means.

read in observation data

contents = np.genfromtxt(data_filename, skip_header = 1, delimiter=’ \n’, dtype=np.str)

n_sources = int(contents[2].split()[0]) # it is at the 3rd line fist element
n_topo = int(contents[579].split()[0]) # the 579th line 1st element

x_locations = np.zeros(n_sources) # zero vector
a_spacings = np.zeros(n_sources)
apparent_resistivity_values = np.zeros(n_sources)
a_locations = np.zeros(n_sources)
b_locations = np.zeros(n_sources) # zero vector
m_locations = np.zeros(n_sources)
n_locations = np.zeros(n_sources)
observed_data =np.zeros(n_sources)
Horizontal_xs = np.zeros(n_topo)
Vertical_ys = np.zeros(n_topo)

content_index = 4

loop over sources at certain range in this case, from 5th row to 577th row

for i in range(0,573): # range(4,576)

start by reading in the source info

   content_index = content_index + 1  # read the next line
   x_location, a_spacing, apparent_resistivity_value = contents[content_index].split() # this is a string(read the 5th line values and assign them to parameters)

convert the strings to a int for locations, ‘a’ spacing and float for apparent resisitivty values

   x_locations[i] = int(x_location)  # mid point location for wenner geometry
   a_spacings[i] = int(a_spacing)
   apparent_resistivity_values[i] = float(apparent_resistivity_value)
   


   a_electrodes = x_locations - (a_spacings/2 + a_spacings)
   b_electrodes = x_locations + (a_spacings/2 + a_spacings)

   m_electrodes = x_locations - (a_spacings/2)
   n_electrodes = x_locations + (a_spacings/2)
   d_obs = apparent_resistivity_values

convert to the UBC format of wenner

a_electrodes = np.vstack([a_electrodes, np.zeros_like(a_electrodes), np.zeros_like(a_electrodes)]).T
b_electrodes = np.vstack([b_electrodes, np.zeros_like(b_electrodes), np.zeros_like(b_electrodes)]).T

m_electrodes = np.vstack([m_electrodes, np.zeros_like(m_electrodes), np.zeros_like(m_electrodes)]).T
n_electrodes = np.vstack([n_electrodes, np.zeros_like(n_electrodes), np.zeros_like(n_electrodes)]).T

Define survey

unique_tx, k = np.unique(np.c_[a_electrodes, b_electrodes], axis=0, return_index=True)
k= np.r_[k,len(a_electrodes)+1] # k = 572, len(a_electrodes) = 573

source_list = []

for ii in range(0,n_sources):

    # MN electrode locations for receivers. Each is an (N, 3) numpy array
    M_locations = m_electrodes[k[ii]:k[ii+1],:]
    N_locations = n_electrodes[k[ii]:k[ii+1],:]
    receiver_list = [dc.receivers.Dipole(M_locations, N_locations)]

    # AB electrode locations for source. Each is a (1, 3) numpy array
    A_location = a_electrodes[k[ii],:]
    B_location = b_electrodes[k[ii],:]
    source_list.append(dc.sources.Dipole(receiver_list, A_location, B_location))

survey = dc.survey.Survey_ky(source_list)

Define the a data object. Uncertainties are added later

dc_data = Data(survey=survey,dobs= d_obs)

plot psuedosection

mpl.rcParams.update({“font.size”: 12})
fig = plt.figure(figsize=(12, 5))

ax1 = fig.add_axes([0.05, 0.05, 0.8, 0.9])
dc.utils.plot_pseudoSection(
dc_data,
ax=ax1,
survey_type=“wenner”,
data_type=“appResistivity”,
space_type=“half-space”,
scale=“log”,
pcolorOpts={“cmap”: “viridis”},
)
ax1.set_title(“Apparent Resistivity [S/m]”)

plt.show()

Thanks!

Hi @Yun,
Sorry for the late reply.
Indeed, it would be easier if you could share:

  1. The data file (you can anonymize the observations if this is sensitive)
  2. your complete script in a .py file

A link to downloadable versions of them could work.

Hi, Thibaut, Thanks, I’ve shared my files to your Email address. [personal information masked.]

Yun

received. Will have a look and let you know.

Hi Thibaut,

My question is very simple and easy to answer,

Learning From the sample code you provided on the public webpage, (2.5D DC Resistivity Least-Square Inversion and 1 century DCIP inversion), I tried to practice doing 2D inversion with python and Simpeg code, I thought it would be easy to do because I can simply follow the instruction you provided to produce any result I want from slightly modifying the code. But when comes to plot the Peseudosection, no matter how similar I set the data to replicate the format of the sample data (even exactly the same, only difference is my data is Wenner array), I couldn’t get it to work, I checked several times, it is particularly dc.utils.plot_pseudosection(), and the problem is how to arrange dc_data as input for it.

Because lack of examples of how to do the 2D Wenner array in Simpeg, I couldn’t get a clear picture of how to do it.
If you don’t have time for this, you can just forward me a simple workable example of how to do a 2D Wenner array, it will give a lot of help for how to arrange the data the right way.

Thanks.
Yun

Ah yes, I see. This is the same issue as the one described in Topo with DC-2d-inversion-app - #5 by thibaut.astic

Short explanation: there is a typo in the Pseudosections plotting function. It is being fixed soon. This bug does not affect at all the inversion (just the plotting of the pseudo-section).

Dis you install SimPEG via pip or conda, or are you working from a github clone? If the latter, you could switch to the dcip_update_utils branch while we update the main branch and the pip and conda install.

Hi
I think it will be useful to creat inputting module to read different type data such as Syscal, ABEM-Lund , E4D, BERT for doing inversing.
Best wishes.

Hi Thibaut, Thanks, I’m still a new user of Simpeg, I installed it from Anaconda which basically covered all the packages I want to work with. I usually keep tracking the update from conda-forge.

Only the plotting of the pseudosection is affected. So feel confident in continuing your project, generating and inverting data. The plotting of the inverted model will have no issue either.

The fix is being brought with a bunch of other various improvements for DC, so that is why it’s taking a little longer. Keep an eye open for new SimPEG releases on Conda, it is indeed a good way to install it :+1:t2: