Creating 2D forward modelling with simpeg saving the outpout plot as matrix or how to save the data of restivity variation wrt x and z axis

DC and IP
1

I want to save the plot of forward model and initial defined model in the matrix of shape N*M for the purpose of training a CNN using the U-NET Machine learning model can anyone help me regarding this…

Also, how can we save the data of resistivity of forward model and initial model with x and z dimensions if this is also possible …

Below is the code for forward model and true model

Thanks,

Code

from SimPEG.electromagnetics.static import resistivity as DC, utils as DCutils
import discretize
from SimPEG import (
maps,
utils,
data_misfit,
regularization,
optimization,
inversion,
inverse_problem,
directives,
)
import matplotlib.pyplot as plt
from matplotlib import colors
import numpy as np
from pylab import hist

try:
from pymatsolver import PardisoSolver as Solver
except ImportError:
from SimPEG import SolverLU as Solver

def run(
plotIt=True,
survey_type=“dipole-dipole”,
rho_background=1e3,
rho_block=1e2,
block_x0=100,
block_dx=10,
block_y0=-10,
block_dy=5,
):

np.random.seed(1)
# Initiate I/O class for DC
IO = DC.IO()
# Obtain ABMN locations

xmin, xmax = 0.0, 200.0
ymin, ymax = 0.0, 0.0
zmin, zmax = 0, 0
endl = np.array([[xmin, ymin, zmin], [xmax, ymax, zmax]])
# Generate DC survey object
survey = DCutils.gen_DCIPsurvey(
    endl, survey_type=survey_type, dim=2, a=10, b=10, n=10
)
survey = IO.from_ambn_locations_to_survey(
    survey.locations_a,
    survey.locations_b,
    survey.locations_m,
    survey.locations_n,
    survey_type,
    data_dc_type="volt",
)

# Obtain 2D TensorMesh
mesh, actind = IO.set_mesh()
# Flat topography
actind = utils.surface2ind_topo(mesh, np.c_[mesh.vectorCCx, mesh.vectorCCx * 0.0])
survey.drape_electrodes_on_topography(mesh, actind, option="top")
# Use Exponential Map: m = log(rho)
actmap = maps.InjectActiveCells(mesh, indActive=actind, valInactive=np.log(1e8))
parametric_block = maps.ParametricBlock(mesh, slopeFact=1e2)
mapping = maps.ExpMap(mesh) * parametric_block
# Set true model
# val_background,val_block, block_x0, block_dx, block_y0, block_dy
mtrue = np.r_[np.log(1e3), np.log(10), 100, 10, -20, 10]

# Set initial model
m0 = np.r_[
    np.log(rho_background),
    np.log(rho_block),
    block_x0,
    block_dx,
    block_y0,
    block_dy,
]
rho = mapping * mtrue
rho0 = mapping * m0
# Show the true conductivity model
fig = plt.figure(figsize=(12, 3))
ax = plt.subplot(111)
temp = rho.copy()
temp[~actind] = np.nan
out = mesh.plotImage(
    temp,
    grid=False,
    ax=ax,
    gridOpts={"alpha": 0.2},
    clim=(10, 1000),
    pcolorOpts={"cmap": "viridis", "norm": colors.LogNorm()},
)
ax.plot(survey.electrode_locations[:, 0], survey.electrode_locations[:, 1], "k.")
ax.set_xlim(IO.grids[:, 0].min(), IO.grids[:, 0].max())
ax.set_ylim(-IO.grids[:, 1].max(), IO.grids[:, 1].min())
cb = plt.colorbar(out[0])
cb.set_label("Resistivity (ohm-m)")
ax.set_aspect("equal")
ax.set_title("True resistivity model")
plt.show()
# Show the true conductivity model
fig = plt.figure(figsize=(12, 3))
ax = plt.subplot(111)
temp = rho0.copy()
temp[~actind] = np.nan
out = mesh.plotImage(
    temp,
    grid=False,
    ax=ax,
    gridOpts={"alpha": 0.2},
    clim=(10, 1000),
    pcolorOpts={"cmap": "viridis", "norm": colors.LogNorm()},
)
ax.plot(survey.electrode_locations[:, 0], survey.electrode_locations[:, 1], "k.")
ax.set_xlim(IO.grids[:, 0].min(), IO.grids[:, 0].max())
ax.set_ylim(-IO.grids[:, 1].max(), IO.grids[:, 1].min())
cb = plt.colorbar(out[0])
cb.set_label("Resistivity (ohm-m)")
ax.set_aspect("equal")
ax.set_title("Initial resistivity model")
plt.show()

# Generate 2.5D DC problem
# "N" means potential is defined at nodes
prb = DC.Simulation2DNodal(
    mesh, survey=survey, rhoMap=mapping, storeJ=True, solver=Solver
)

# Make synthetic DC data with 5% Gaussian noise
data = prb.make_synthetic_data(mtrue, relative_error=0.05, add_noise=True)

# Show apparent resisitivty pseudo-section
IO.plotPseudoSection(data=data.dobs / IO.G, data_type="apparent_resistivity")
##

if name == “main”:
run()
plt.show()

2

Hi @priyamomer ,

The model is a vector ordered in the same way as mesh.gridCC
so you can access the coordinates through mesh.gridCC and the model values through the model.
After that, these are all numpy object so can use the reshape function to order it the way you want.