# DC-2d inversion with Wenner

**URL:** <https://simpeg.discourse.group/t/dc-2d-inversion-with-wenner/196>\
**Category:** DC and IP\
**Created:** [March 15, 2021, 7:51am UTC](https://simpeg.discourse.group/t/dc-2d-inversion-with-wenner/196 "2021-03-15T07:51:01Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![Yun](https://avatars.discourse-cdn.com/v4/letter/y/bcef8e/32.png) [@Yun](https://simpeg.discourse.group/u/Yun)\
**Post date:** [March 15, 2021, 7:51am UTC](https://simpeg.discourse.group/t/dc-2d-inversion-with-wenner/196/1 "2021-03-15T07:51:01Z")

</div>

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](https://docs.simpeg.xyz/content/tutorials/05-dcr/plot_inv_2_dcr2d.html) 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.

---

<div class="post-metadata">

**Author:** ![thibaut.astic](https://yyz2.discourse-cdn.com/free1/user_avatar/simpeg.discourse.group/thibaut.astic/32/109_2.png) [@thibaut.astic](https://simpeg.discourse.group/u/thibaut.astic)\
**Post date:** [March 17, 2021, 5:07am UTC](https://simpeg.discourse.group/t/dc-2d-inversion-with-wenner/196/2 "2021-03-17T05:07:24Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Yun](https://avatars.discourse-cdn.com/v4/letter/y/bcef8e/32.png) [@Yun](https://simpeg.discourse.group/u/Yun)\
**Post date:** [March 18, 2021, 12:50am UTC](https://simpeg.discourse.group/t/dc-2d-inversion-with-wenner/196/3 "2021-03-18T00:50:00Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![thibaut.astic](https://yyz2.discourse-cdn.com/free1/user_avatar/simpeg.discourse.group/thibaut.astic/32/109_2.png) [@thibaut.astic](https://simpeg.discourse.group/u/thibaut.astic)\
**Post date:** [March 26, 2021, 5:24pm UTC](https://simpeg.discourse.group/t/dc-2d-inversion-with-wenner/196/4 "2021-03-26T17:24:27Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Yun](https://avatars.discourse-cdn.com/v4/letter/y/bcef8e/32.png) [@Yun](https://simpeg.discourse.group/u/Yun)\
**Post date:** [April 5, 2021, 4:26pm UTC](https://simpeg.discourse.group/t/dc-2d-inversion-with-wenner/196/5 "2021-04-05T16:26:17Z")

</div>

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

Yun

---

<div class="post-metadata">

**Author:** ![thibaut.astic](https://yyz2.discourse-cdn.com/free1/user_avatar/simpeg.discourse.group/thibaut.astic/32/109_2.png) [@thibaut.astic](https://simpeg.discourse.group/u/thibaut.astic)\
**Post date:** [April 5, 2021, 8:03pm UTC](https://simpeg.discourse.group/t/dc-2d-inversion-with-wenner/196/6 "2021-04-05T20:03:03Z")

</div>

received. Will have a look and let you know.

---

<div class="post-metadata">

**Author:** ![Yun](https://avatars.discourse-cdn.com/v4/letter/y/bcef8e/32.png) [@Yun](https://simpeg.discourse.group/u/Yun)\
**Post date:** [April 6, 2021, 11:24pm UTC](https://simpeg.discourse.group/t/dc-2d-inversion-with-wenner/196/7 "2021-04-06T23:24:20Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![thibaut.astic](https://yyz2.discourse-cdn.com/free1/user_avatar/simpeg.discourse.group/thibaut.astic/32/109_2.png) [@thibaut.astic](https://simpeg.discourse.group/u/thibaut.astic)\
**Post date:** [April 14, 2021, 10:06pm UTC](https://simpeg.discourse.group/t/dc-2d-inversion-with-wenner/196/8 "2021-04-14T22:06:28Z")

</div>

Ah yes, I see. This is the same issue as the one described in [Topo with DC-2d-inversion-app - #5 by thibaut.astic](https://simpeg.discourse.group/t/topo-with-dc-2d-inversion-app/183/5)

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.

---

<div class="post-metadata">

**Author:** ![subenyu](https://avatars.discourse-cdn.com/v4/letter/s/e19b73/32.png) [@subenyu](https://simpeg.discourse.group/u/subenyu)\
**Post date:** [April 17, 2021, 10:16pm UTC](https://simpeg.discourse.group/t/dc-2d-inversion-with-wenner/196/9 "2021-04-17T22:16:45Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Yun](https://avatars.discourse-cdn.com/v4/letter/y/bcef8e/32.png) [@Yun](https://simpeg.discourse.group/u/Yun)\
**Post date:** [April 19, 2021, 6:56am UTC](https://simpeg.discourse.group/t/dc-2d-inversion-with-wenner/196/10 "2021-04-19T06:56:17Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![thibaut.astic](https://yyz2.discourse-cdn.com/free1/user_avatar/simpeg.discourse.group/thibaut.astic/32/109_2.png) [@thibaut.astic](https://simpeg.discourse.group/u/thibaut.astic)\
**Post date:** [April 23, 2021, 5:23pm UTC](https://simpeg.discourse.group/t/dc-2d-inversion-with-wenner/196/11 "2021-04-23T17:23:08Z")

</div>

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 👍🏻

---

<div class="post-metadata">

**Author:** ![vrath](https://yyz2.discourse-cdn.com/free1/user_avatar/simpeg.discourse.group/vrath/32/206_2.png) [@vrath](https://simpeg.discourse.group/u/vrath)\
**Post date:** [September 5, 2022, 4:03pm UTC](https://simpeg.discourse.group/t/dc-2d-inversion-with-wenner/196/12 "2022-09-05T16:03:55Z")

</div>

Hi,  
I totally agree! Maybe somebody in the community has experience with Syscal .bin files, and can share a module/routine importing?  
Thanks, Volker
