Package Modeling

Module Occam2D

This file is part of pyCSAMT.

pyCSAMT is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

pyCSAMT is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License along with pyCSAMT. If not, see <https://www.gnu.org/licenses/>.

---

See also

• http://marineemlab.ucsd.edu/Projects/Occam/sharp/index.html

• https://marineemlab.ucsd.edu/~kkey/Software.php

  …

Note

Actually this module is not able to write Occam2D meshes file or build OccamInputfiles. However , you can import ocam module of :ref: MTpy toolbox …

@author: KouaoLaurent alias @Daniel03 Created on Fri Jan 22 20:31:14 2021

class pycsamt.modeling.occam2d.Data(data_fn=None, **kwargs)

Read and write Occam data

Methods

parse_block_data([datablocks])	From block data, retreive the value of each blocks each blocks is ndaray(nfrea, nstations) Attributes are : - self.appRho for apparent resistivity blocks - self.phase for phase values - self.error_appRho for error in resistivity (block ) - self.error_phase for error in phase 9blocks
read_occam2d_datafile([data_fn])	read_occam_data file.

parse_block_data(datablocks=None)

From block data, retreive the value of each blocks each blocks is ndaray(nfrea, nstations) Attributes are :

• self.appRho for apparent resistivity blocks

• self.phase for phase values

• self.error_appRho for error in resistivity (block )

• self.error_phase for error in phase 9blocks

read_occam2d_datafile(data_fn=None)

read_occam_data file. and populates attributes

Parameters

fn (str) – full path to occam data file .

Example

>>> path =os.path.join(os.environ ['pyCSAMT'],
...                       'pycsamt', 'data', 'occam2D',
...                       'OccamDataFile.dat')
>>> occamdata_obj =Data(fn = path)
>>> print(occamdata_obj.occam_data_format)
>>> print(occamdata_obj.occam_data_blocks)
>>> print(occamdata_obj.occam_data_title)
>>> print(occamdata_obj.occam_data_sites)
>>> print(occamdata_obj.occam_data_frequencies)
>>> print(occamdata_obj.occam_data_offsets)
>>> occamDATA = occamdata_obj.occam_data

class pycsamt.modeling.occam2d.DataBlock(**kwargs)

Read OccamDataBlock aand set corresponding attributes

Attributes	Description (Restriction)
site	number of the site from the site list that this data belongs to.
freq	number of the frequency from the frequency list.
type	type of data
datum	data value
error	size of the error for this measurement For log10 resistivity this value can look a little strange. It is derived from the calculation d(log_10(x))/dx = 1/[x ln(10)]. So for 10% error, dx = 0.10x thus d(log_10(x)) = 0.10x / x ln(10) =.1/ln(10) = 0.0434

Example

>>> form pycsamt.modeling import DataBlock
>>> np.random.seed(1983)
>>> dblocks_dict = {'site' : ['S{0:02}'.format(ii)  for ii in range (10)],
...                'freq': np.linspace(1, 8912, 10),
...                'type' :(5,6),
...                'data' : np.power(10, np.random.randn(20*2)),
...                'error' : np.random.randn(20)
...                }
>>> dblock_obj = DataBlock(**dblocks_dict)
>>> print(dblock_obj.freq )
>>> print(dblock_obj.error)
>>> print(dblock_obj.site)

Methods

decode_each_site_data(data_blocks, ...)

Decode each site data and get differents values array :param datablocks: datafrom differents blocks :type datablocks: ndarray

static decode_each_site_data(data_blocks, data_type_index)

Decode each site data and get differents values array :param datablocks: datafrom differents blocks :type datablocks: ndarray

Parameters

data_type_index (int) – specify the data type index

class pycsamt.modeling.occam2d.Iter(iter_fn=None, **kwargs)

Occam iteration file, inherets from startup obj , in fact two object a similar the same

Methods

read_occam2d_iterfile([iter_fn])	read_occam iteration file and populate attributes
read_occam2d_startupfile([startup_fn])	read occam2d_startupfile

read_occam2d_iterfile(iter_fn=None)

read_occam iteration file and populate attributes

Parameters

iter_fn (str) – full path to iteration file

Example

>>> from pycsamt.modeling.occam2d import Iter
>>> path =os.path.join(os.environ ['pyCSAMT'], 'pycsamt',
...                       'data', 'occam2D', 'ITER17.iter')
>>> iter_obj=Iter(iter_fn= path)
... print(iter_obj.occam_iter_param_count)
... iterDATA= iter_obj.occam_iter_data

class pycsamt.modeling.occam2d.Iter2Dat(iter2dat_fn=None, model_fn=None, data_fn=None, iter_fn=None, mesh_fn=None, **kwargs)

Iter2Dat is a format converter which convert *.iter file and related mesh files to so called ‘x,y,z’ *.data file for post-processing. The class was inpired from the Bo Yang matlab script . reading functions come from ‘plotOccam2DMT.m’ routine

See also

Occam routine <http://marineemlab.ucsd.edu/Projects/Occam/sharp/index.html > of SIO, UCSD.

Bo Yang matlab-script is on add.info sub-packages of pyCSAMT. If your are familiar with matlab and you want to rewrite the script please contact author at:

Methods

read_iter2dat([iter2dat_fn, bln_fn, scale, ...])	Read YangBo iter2data file or provided Occam 2D specific files .
read_occamfiles([path_to_occamfiles])	getffiles and readfiles to populates speciales attributes
write_iter2dat_file([filename, model_fn, ...])	write 'x,y,z' *.data file for post-processing can read and rewrite iter2dat file

read_iter2dat(iter2dat_fn=None, bln_fn=None, scale='km', model_fn=None, iter_fn=None, mesh_fn=None, doi='1km', data_fn=None, occam_model_obj=None, **kws)

Read YangBo iter2data file or provided Occam 2D specific files .

Parameters

• iter2dat_fn (str) – full path to iter2dat file

• bln_file (str) – full path to bln file

• scale –

  str , scale of output data . Most of time , Bo yang iter2Dat file is converted in kilometer . If not turn

  the scale to m.

Example

>>> from pycsamt.modeling.occam2d import Iter2Dat
>>> i2d='iter17.2412.dat'
>>> i2d_2='K1.iter.dat'
>>> bln='iter17.2412.bln'
>>> bln_2 = 'K1.bln'
>>> pathiter = os.path.join(os.path.dirname(os.environ ['pyCSAMT']),
...                            '_iter2dat_', i2d_2)
>>> pathbln = os.path.join(os.path.dirname(os.environ ['pyCSAMT']),
...                           '_iter2dat_', bln_2)
>>> occam_iter2dat_obj =Iter2Dat(iter2dat_fn=pathiter,
...                                 bln_fn =pathbln )
>>> i2d_data = occam_iter2dat_obj.read_iter2datfile()
>>> i2d_sta, i2d_depth = occam_iter2dat_obj.model_x_nodes,
>>> occam_iter2dat_obj.model_z_nodes
>>> i2d_model_res = occam_iter2dat_obj.model_res

read_occamfiles(path_to_occamfiles=None, **kws)

getffiles and readfiles to populates speciales attributes

Parameters

path_to_occamfiles (str) – full path to occamfiles [ITER|DATA|RESP] files.

Example

>>> path =os.path.join(os.environ ['pyCSAMT'], 'pycsamt', 'data', 'occam2D')
>>> iter2_obj=Iter2Dat(path_to_occamfiles= path)
>>> sites= iter2_obj.OccamData.occam_data_sites
>>> freq =iter2_obj.OccamData.occam_data_frequencies
>>> iter_roughn =iter2_obj.OccamIter.occam_iter_roughness_value
>>> iter_misfit=iter2_obj.OccamIter.occam_iter_misfit_reached
>>> forward_data = iter2_obj.OccamResponse.forward_data
>>> residual_data = iter2_obj.OccamResponse.residual

write_iter2dat_file(filename=None, model_fn=None, iter_fn=None, mesh_fn=None, data_fn=None, doi='1km', savepath=None, occam_model_obj=None, negative_depth=True, scale='km', **kws)

write ‘x,y,z’ *.data file for post-processing can read and rewrite iter2dat file

Parameters

• x (array_like) – offset

• y (array_like) – depth

• z (array_like) – log10 resistivities

params	Type	Description
model_fn	str	full path to occam model file
iter_fn	str	full path to iteration file
mesh_fn	str	full path to occam mesh file
data_fn	str	fll path to occam data file
filename	str	output of iter2dat file
doi	str of float	depth of investigation. default is 1km if your povide value on float, default unit is “meter” eg : 2000=2000m
negative_depth	bool	if True , will provide file with negative depth value
scale	str	scaled the offset value and elevation . might be [m|km]
elevation	ndarray|list	can be provided if usefull

1. Write with Occam 2D files

Examples

>>> from pycsamt.modeling.occam2d import Iter2Dat as i2d
>>> data='OccamDataFile.dat'
>>> mesh = 'Occam2DMesh'
>>> model = 'Occam2DModel'
>>> path =os.path.join(os.environ ['pyCSAMT'], 'pycsamt',
...                       'data', 'occam2D')
...                   #,'OccamDataFile.dat')
>>> pathi2d =os.path.join(os.environ ['pyCSAMT'], 'pycsamt', 'data', '_iter2dat_')
>>> occam_iter2dat_obj =i2d(mesh_fn=os.path.join(path, mesh),
...                    iter_fn = os.path.join(path, iter_),
...                    model_fn =os.path.join(path, model),
...                    data_fn =os.path.join(path, data))
>>> occam_iter2dat_obj.write_iter2dat_file()

2. Rewrite the with iter2dat file

Example

>>> from pycsamt.modeling.occam2d import Iter2Dat as i2d
>>> iter_='ITER17.iter'
>>> idat = K1.iter.20142.dat
>>> bln = K1.iter.20142.bln
>>> occam_iter2dat_obj =i2d(iter2dat_fn=os.path.join(pathi2d, idat),
                            bln_fn = os.path.join(pathi2d, bln)
>>> occam_iter2dat_obj.write_iter2dat_file()

class pycsamt.modeling.occam2d.Mesh(mesh_fn=None, **kwargs)

Read Occam read mesh file

Methods

read_occam2d_mesh([mesh_fn])

Read mesh and get mesh values data and populates attributes

read_occam2d_mesh(mesh_fn=None)

Read mesh and get mesh values data and populates attributes

Parameters

mesh_fn (str) – full path to mesh file

class pycsamt.modeling.occam2d.Model(model_fn=None, iter_fn=None, mesh_fn=None, **kwargs)

Occam Model , Actually read model file

Methods

read_occam2d_modelfile([model_fn, mesh_fn, ...])

read and ascertain modelfile.

read_occam2d_modelfile(model_fn=None, mesh_fn=None, iter_fn=None)

read and ascertain modelfile.

Parameters

• model_fn (str) – full path to OCCAM2D model file

• mesh_fn (str) – full path to MESH model file

• iter_fn (str) – full path to ITER model file

Example

>>> path =os.path.join(os.environ ['pyCSAMT'], 'pycsamt',
...                       'data', 'occam2D', 'Occam2DModel')
>>> occammodel_obj =Model(model_fn = path)
... print(occammodel_obj.model_name)
... print(occammodel_obj.model_mesh_file)
... print(occammodel_obj.model_num_layers)

class pycsamt.modeling.occam2d.Response(response_fn=None, data_fn=None, **kwargs)

Occam response file . Response is paired with the iteration file above. One is output at the end of each successful iteration. Inherets from Occam 2D data file .

Methods

parse_block_data([datablocks])	From block data, retreive the value of each blocks each blocks is ndaray(nfrea, nstations) Attributes are : - self.appRho for apparent resistivity blocks - self.phase for phase values - self.error_appRho for error in resistivity (block ) - self.error_phase for error in phase 9blocks
read_occam2d_datafile([data_fn])	read_occam_data file.
read_occam2d_responsefile([response_fn])	Read Occam2D response file.

read_occam2d_responsefile(response_fn=None)

Read Occam2D response file.

Parameters

response_fn – full path to response file

:type response_fn:str

Example

>>> from pycsamt.modeling.occam2d import Response
>>> pathresp =os.path.join(os.environ ['pyCSAMT'], 'pycsamt',
...                           'data', 'occam2D',RESP17.resp )
>>> path_data =os.path.join(os.environ ['pyCSAMT'], 'pycsamt',
...                            'data', 'occam2D',OccamDataFile.dat )
>>> resp_obj = Response(response_fn=pathresp, data_fn = path_data )
>>> respDATA= resp_obj.resp_data_value
>>> resp_obj.occam_mode
>>> resp_obj.occam_dtype
>>> forward_data = resp_obj.resp_forward_value
>>> residual_data = resp_obj.resp_residual_value
>>> tm_log10= resp_obj.resp_tm_log10
>>> tm_phase=resp_obj.resp_tm_phase
>>> tm_forward = resp_obj.resp_tm_log10_forward
>>> tm_residual = resp_obj.resp_tm_log10_forward
>>> tm_forward = resp_obj.resp_tm_log10_residual
>>> tm_phase_error = resp_obj.resp_tm_phase_err
>>> tm_phase_err = resp_obj.resp_tm_log10_forward_phase
>>> tm_residual_phase_err = resp_obj.resp_tm_log10_residual_phase_err

class pycsamt.modeling.occam2d.Startup(startup_fn=None, **kwargs)

Occam startup file Actually read startup file. for more detail:

See also

http://marineemlab.ucsd.edu/Projects/Occam/sharp/index.html

Methods

read_occam2d_startupfile([startup_fn])

read occam2d_startupfile

read_occam2d_startupfile(startup_fn=None)

read occam2d_startupfile

Parameters

startup_fn (str) – full path to startup_file

Example

>>> from pycsamt.modeling.occam2d import Startup
>>> path =os.path.join(os.environ ['pyCSAMT'],
...                       'pycsamt', 'data', 'occam2D', 'Startup')
>>> startup_obj=Startup(startup_fn = path)
... print(startup_obj.occam_startup_data_file)

pycsamt.modeling.occam2d.getMisfit(resp_fn=None, data_fn=None, kind='rho', **kwargs)

Calling getMisfit to plot misfit value using geoplot2d decorator see pycsamt.viewer.plot.geoplot2d to get the kwargs arguments. Most keywords arguments are the same used by pycsamt.viewer.plot.Plot2d Please refer to mod documentation. For more details , refer to :class:pycsamt.viewer.plot.geoplot2d

Warning:

when using decorator`geoplot2d`, set always reason argument to misfit or else.

Parameters

kwargs (dict) – use the main argument from getMisfit function

:return:* resp_misfit,

• resp_sites_names,

• resp_sites_offsets,

• resp_freq, doi ,

• model_rms,

• model_roughness

Example

>>> from pycsamt.modeling.occamd2d import getMisfit
>>> occamPath ='data/occam2d
>>> response_file ='RESP27.resp'
>>> data_file = 'OccamDataFile.dat'
>>> log_file = 'LogFile.logfile'
>>> getMisfit(response_fn = os.path.join(occamPath, response_file)
...            data_fn = os.path.join(occamPath, data_file),
...            logfile =  os.path.join(occamPath, log_file))

class pycsamt.modeling.occam2d.occam2d_write

Special class to build occam2d imput files with MTpy module .

Parameters

• **edi_fn** (str) – full path to edifiles locations

• **freq_num** (float) – number of frequencies to use in inversion

• **interpolate_freq** (bool,) – frequency interpolation , default is False

• **geoelectric_strike** (bool) – geoelectric strike angle assuming N = 0, E = 90. If True , provided , losgspace interpolation as tuple value

Others important attributes can be found in :

Key Words/Attributes	Description
edi_fn	full path to data file
n_layers	number of vertical layers in mesh default is 31.
num_layers	[ int ] number of regularization layers.
num_z_pad_cells	number of vertical padding cells below
iterations_to_run	maximum number of iterations to run default is 20
resistivity_start	starting resistivity value. If model_values is not given, then all values with in model_values array will be set to resistivity_start
save_path	directory path to save startup file to default is current working directory
startup_basename	basename of startup file name. default is Occam2DStartup
startup_fn	full path to startup file. default is save_path/startup_basename
target_misfit	target misfit value. default is 1.
x_pad_multiplier	horizontal padding cells will increase by this multiple out to the edge of the grid. default is 1.7
z1_layer	thickness of the first layer in the model. Should be at least 1/4 of the first skin depth default is 10
z_bottom	bottom depth of the model (m). Needs to be large enough to be 1D at the edge. default is 200000.0
z_target_depth	depth to deepest target of interest. Below this depth cells will be padded to z_bottom
cell_width	width of cells with in station area in meters default is 100
phase_te_err	percent error in phase for TE mode. default is 5
phase_tm_err	percent error in phase for TM mode. default is 20.
res_te_err	percent error in resistivity for TE mode. default is 10
res_tm_err	percent error in resistivity for TM mode.
trigger	[ float ] multiplier to merge model blocks at depth. A higher number increases the number of model blocks at depth. default is .1.12
model_mode	model mode to use for inversion, see module`Data`.

Note

We consider occam2D buildingInputs file is focused on CSAMT data in TM mode then default configuration as setting according this feature. If input edi_fn are MT data , please resetting configuration you SEG edi-fn data provided !. …

Example

>>> from pycsamt.modeling.occam2d import occam2d_write
>>>  occam2d_write.buildingInputfiles(os.path.join(os.environ['pyCSAMT'],
...                                              'data', 'edi'),
...                        savepath =os.path.join(os.environ['pyCSAMT'],
...                                               'data', 'tesocc2' ),
...                        geolectricke_strike=34.)

Methods

buildingInputfiles(edi_fn[, freq_num, ...])

Method to build Occam2D inputfiles.

static buildingInputfiles(edi_fn, freq_num=None, savepath=None, interpolate_freq=False, geoelectric_strike=None, **kwargs)

Method to build Occam2D inputfiles. Deal with MTpy module. Try to install MTpy is not installed yet.

Parameters

• edi_fn (str) – full path to edifiles

• freq_num (float) – number of frequencies to use in inversion

• interpolate_freq (bool) – frequency interpolation , default is False

• geoelectric_strike (float) – geoelectric strike angle assuming N = 0, E = 90. If True , provided , losgspace interpolation as tuple value

Returns

outfiles building occam2d inputfiles

Return type

str , sys.stdout

class pycsamt.modeling.occam2d.occamLog(fn=None, **kwargs)

Class to deal with occcam 2d logfile . output File after inverted data. most likely called logFile.logfile.

Attributes

fn

Methods

read_occam2d_logfile([fn])

Read occam file and populate attributes.

read_occam2d_logfile(fn=None)

Read occam file and populate attributes.

Parameters

fn (str) – full path to occam2d log file

Example

>>> from pycsamt.modeling import occam2d
>>> path =os.path.join(os.environ ['pyCSAMT'], 'pycsamt', 'data',
...                       'occam2D', 'logFile.logfile')
>>> occamlog_obj =occam2d.Log(fn = path)
>>> print(occamlog_obj.rms.shape)
>>> print(occamlog_obj.roughness.shape)

pycsamt.modeling.occam2d.plotResponse(data_fn=None, resp_fn=None, stations=None, **kws)

Find Responses properties and plot values and could read multiples lines provided that each line as its occam datafile and its response file.

Response inherits of Data class then easy to data attributes for plot purposes. Visualize all error from raw to inversion using error_type arguments. For instance:

- `1` : visualize the misfit compute manually
- `2`:  visualize the raw error from raw occam data
- `3`: visualzie the error data and phase  defined as
        * error = (input data - foward data)/ RESI
- `4` or 'residual`': Visualize only  the residual data .

Default is `residual.`

Parameters

• data_fn – Occam data file. Can be a string or list of datafiles

• resp_fn – Response file and can be string or list of response file.

• stations – str or list of station to visualize.

Note

To visualize multiple files. Can prodile only the path

where occam (data and response) files are located. In that case, set datafile and response file the name like:

`data_fn` =['line01.dat' , 'line01.resp', 'line02.dat',
          'line02.resp' , ...]

Returns

• resp_lines : coloection of survey lines id

• resp_stations : Collection of station to visualize

• resp_freq: collection of frequencies

• resp_appRHO: collection of te/tm data and

  forward data as a tuple(te|tm , fw’te|tm’)

• resp_phase: collection of the te/tm phase data as tuple (data, fw)

• resp_appRho_err. collection of misfit rho response

• resp_phase_err: collection of the misfit phase

Example

>>> from pysamt.modeling.occam2d import plotResponse
>>> resPath =r'F:\ThesisImp\occam2D\old\K1'
>>> plotResponse(data_fn =resPath,
...                 stations = ['S00', 'S04', 's08', 'S12']
...                  rms =['1.013', '1.451', '1.00', '1.069'],
...                  error_type ='resi' )