mtpy.modeling.modem package

Submodules

mtpy.modeling.modem.config module

ModEM

# Generate files for ModEM

# revised by JP 2017 # revised by AK 2017 to bring across functionality from ak branch

class mtpy.modeling.modem.config.ModEMConfig(**kwargs)[source]

Bases: object

Read and write configuration files for how each inversion is run.

add_dict(fn=None, obj=None)[source]: Add dictionary based on file name or object.

write_config_file(save_dir=None, config_fn_basename='ModEM_inv.cfg')[source]: Write a config file based on provided information.

mtpy.modeling.modem.control_fwd module

ModEM

# Generate files for ModEM

# revised by JP 2017 # revised by AK 2017 to bring across functionality from ak branch # revised by OA 2024 to update docs

class mtpy.modeling.modem.control_fwd.ControlFwd(**kwargs)[source]

Bases: object

Reads and writes control files for ModEM to control how the forward solver starts and how it is run.

num_qmr_iter Number of QMR iterations per divergence correction.: int, default is 40
max_num_div_calls Maximum number of divergence correction calls. int,: default is 20
max_num_div_iters Maximum number of divergence correction iterations.: int, default is 100
misfit_tol_fwd Misfit tolerance for EM forward solver. float,: default is 1e-07
misfit_tol_adj Misfit tolerance for EM adjoint solver. float,: default is 1e-07
misfit_tol_div Misfit tolerance for divergence correction. float,: default is 1e-05
save_path Path at which to save the control file. PosixPath,: default is current working directory.
fn_basename Name of the control file. str, default is: ‘control.fwd’,

property control_fn: Control fn.

read_control_file(control_fn=None)[source]: Read in a control file.

write_control_file(control_fn=None, save_path=None, fn_basename=None)[source]

Write control file.

Arguments::

control_fnstring: full path to save control file to default is save_path/fn_basename
save_pathstring: directory path to save control file to default is cwd
fn_basenamestring: basename of control file default is control.inv

mtpy.modeling.modem.control_inv module

ModEM

# Generate files for ModEM

# revised by JP 2017 # revised by AK 2017 to bring across functionality from ak branch # revised by OA 2024 to update docs

class mtpy.modeling.modem.control_inv.ControlInv(**kwargs)[source]

Bases: object

Reads and writes control files for ModEM to control how the inversion starts and how it is run.

output_fn Model and data output file name “prefix phrase”, i.e.: written to the front of ModEM output file names before the iteration number and file extension. str, default is ‘MODULAR_NLCG’

lambda_initial Initial damping factor lambda. int, default is 10 lambda_step To update lambda, divide by this value. int, default

is 10

model_search_step Initial search step in model units. int, default is: 1
rms_reset_search Restart when rms diff is less than this value. float,: default is 0.002
rms_target Exit search when rms is less than this value. float,: default is 1.05
lambda_exit Exit when lambda is less than this value. float,: default is 0.0001

max_iterations Maximum number of iterations. int, default is 100 save_path Path at which to save the control file. PosixPath,

default is current working directory

fn_basename Name of the control file. str, default is: ‘control.inv’

property control_fn: Control fn.

read_control_file(control_fn=None)[source]: Read in a control file.

write_control_file(control_fn=None, save_path=None, fn_basename=None)[source]

Write control file.

Arguments::

control_fnstring: full path to save control file to default is save_path/fn_basename
save_pathstring: directory path to save control file to default is cwd
fn_basenamestring: basename of control file default is control.inv

mtpy.modeling.modem.convariance module

ModEM

# Generate files for ModEM

# revised by JP 2017 # revised by AK 2017 to bring across functionality from ak branch

class mtpy.modeling.modem.convariance.Covariance(grid_dimensions=None, **kwargs)[source]

Bases: object

Class that deals with model covariance and smoothing, writing covariance (.cov) files for use in ModEM.

grid_dimensions Grid dimensions excluding air layers (Nx, Ny, NzEarth) smoothing_east Smoothing in east direction. float, default is 0.3 smoothing_north Smoothing in north direction. float, default is 0.3 smoothing_z Smoothing in vertical direction. float, default is

0.3

smoothing_num Number of smoothing iterations. int, default is 1 exception_list List of exceptions. list, default is empty mask_arr Mask array. numpy.ndarray, default is None save_path Path at which to save the covariance file. PosixPath,

default is current working directory

fn_basename Name of the covariance file. default is: ‘covariance.cov’

property cov_fn: Cov fn.

get_parameters()[source]: Get parameters.

read_cov_file(cov_fn)[source]: Reads a ModEM covariance (.cov) file. :param cov_fn: Filename of the target ModEM covariance (.cov) file. :type cov_fn: str :raises CovarianceError: Error related to the covariance class.

write_cov_vtk_file(cov_vtk_fn, model_fn=None, grid_east=None, grid_north=None, grid_z=None)[source]: Write a vtk file of the covariance to match things up.

write_covariance_file(cov_fn=None, save_path=None, fn_basename=None, res_model=None, sea_water=0.3, air=1000000000000.0)[source]

Writes a ModEM covariance (.cov) file. :param cov_fn: Name for the covariance file. default is None, defaults to None. :type cov_fn: str | Path, optional :param save_path: Location at which to save the covariance file.

default is None, defaults to None.

Parameters:

fn_basename (_type_, optional) – _description_, default is None, defaults to None.
res_model (_type_, optional) – Resistivity model the covariance should be generated from, default is None, defaults to None.
sea_water (float, optional) – Electrical resistivity of sea water cells in Ohm-meters within the model domain, default is 0.3, defaults to 0.3.
air (float, optional) – Electrical resistivity of air cells in Ohm-meters within the model domain, default is 1.0e12, defaults to 1.0e12.

Raises:

CovarianceError – Error related to the covariance class.

mtpy.modeling.modem.data module

ModEM

# Generate files for ModEM

# revised by JP 2017 # revised by AK 2017 to bring across functionality from ak branch # revised by JP 2021 adding functionality and updating. # revised by JP 2022 to work with new structure of a central object

class mtpy.modeling.modem.data.Data(dataframe=None, center_point=None, **kwargs)[source]

Bases: object

Data will read and write .dat files for ModEM and convert a WS data file to ModEM format.

..note: :: the data is interpolated onto the given periods such that all: stations invert for the same periods. The interpolation is a linear interpolation of each of the real and imaginary parts of the impedance tensor and induction tensor. See mtpy.core.mt.MT.interpolate for more details

Parameters:

**kwargs –
center_point – Defaults to None.
dataframe – Defaults to None.
edi_list – List of edi files to read.

property data_filename: Data filename.

property dataframe: Dataframe function.

fix_data_file(fn=None, n=3)[source]: A newer compiled version of Modem outputs slightly different headers This aims to convert that into the older format :param fn: DESCRIPTION, defaults to None. :type fn: TYPE, optional :param n: DESCRIPTION, defaults to 3. :type n: TYPE, optional :return: DESCRIPTION. :rtype: TYPE

get_n_stations(mode)[source]: Get n stations.

property model_parameters: Model parameters.

property n_periods: N periods.

property period: Period function.

read_data_file(data_fn)[source]: Read data file. :param data_fn: Full path to data file name. :type data_fn: string or Path :raises ValueError: If cannot compute component.

write_data_file(file_name=None, save_path=None, fn_basename=None, elevation=False)[source]

Write data file. :param file_name:

Defaults to None.

Parameters:

save_path (string or Path, optional) – Full directory to save file to, defaults to None.
fn_basename (string, optional) – Basename of the saved file, defaults to None.
elevation (boolean, optional) – If True adds in elevation from ‘rel_elev’ column in data array, defaults to False.

Raises:

NotImplementedError – If the inversion mode is not supported.
ValueError – mtpy.utils.exceptions.ValueError if a parameter.

Returns:

Full path to data file.

Return type:

Path

mtpy.modeling.modem.exception module

ModEM

# Generate files for ModEM

# revised by JP 2017 # revised by AK 2017 to bring across functionality from ak branch

exception mtpy.modeling.modem.exception.DataError[source]

Bases: ModEMError

Raise for ModEM Data class specific exceptions.

exception mtpy.modeling.modem.exception.ModEMError[source]: Bases: Exception

mtpy.modeling.modem.model module

ModEM

# Generate files for ModEM

# revised by JP 2017 # revised by AK 2017 to bring across functionality from ak branch # revised by JP 2021 updating functionality and updating docs

class mtpy.modeling.modem.model.Model(station_locations=None, center_point=None, **kwargs)[source]

Bases: object

Make and read a FE mesh grid

The mesh assumes the coordinate system where:: x == North y == East z == + down

All dimensions are in meters.

The mesh is created by first making a regular grid around the station area, then padding cells are added that exponentially increase to the given extensions. Depth cell increase on a log10 scale to the desired depth, then padding cells are added that increase exponentially..

Parameters:: **station_object** – mtpy.modeling.modem.Stations object .. seealso:: mtpy.modeling.modem.Stations

Examples

Example 1 –> create mesh first then data file:

>>> import mtpy.modeling.modem as modem
>>> import os
>>> # 1) make a list of all .edi files that will be inverted for
>>> edi_path = r"/home/EDI_Files"
>>> edi_list = [os.path.join(edi_path, edi)

for edi in os.listdir(edi_path)

>>> ...         if edi.find('.edi') > 0]
>>> # 2) Make a Stations object
>>> stations_obj = modem.Stations()
>>> stations_obj.get_station_locations_from_edi(edi_list)
>>> # 3) make a grid from the stations themselves with 200m cell spacing
>>> mmesh = modem.Model(station_obj)
>>> # change cell sizes
>>> mmesh.cell_size_east = 200,
>>> mmesh.cell_size_north = 200
>>> mmesh.ns_ext = 300000 # north-south extension
>>> mmesh.ew_ext = 200000 # east-west extension of model
>>> mmesh.make_mesh()
>>> # check to see if the mesh is what you think it should be
>>> msmesh.plot_mesh()
>>> # all is good write the mesh file
>>> msmesh.write_model_file(save_path=r"/home/modem/Inv1")
>>> # create data file
>>> md = modem.Data(edi_list, station_locations=mmesh.station_locations)
>>> md.write_data_file(save_path=r"/home/modem/Inv1")

Example 2 –> Rotate Mesh:

>>> mmesh.mesh_rotation_angle = 60
>>> mmesh.make_mesh()

Note

ModEM assumes all coordinates are relative to North and East, and does not accommodate mesh rotations, therefore, here the rotation is of the stations, which essentially does the same thing. You will need to rotate you data to align with the ‘new’ coordinate system.

Attributes	Description
_logger	python logging object that put messages in logging format defined in logging configure file, see MtPyLog more information
cell_number_ew	optional for user to specify the total number of sells on the east-west direction. default is None
cell_number_ns	optional for user to specify the total number of sells on the north-south direction. default is None
cell_size_east	mesh block width in east direction default is 500
cell_size_north	mesh block width in north direction default is 500
grid_center	center of the mesh grid
grid_east	overall distance of grid nodes in east direction
grid_north	overall distance of grid nodes in north direction
grid_z	overall distance of grid nodes in z direction
model_fn	full path to initial file name
model_fn_basename	default name for the model file name
n_air_layers	number of air layers in the model. default is 0
n_layers	total number of vertical layers in model
nodes_east	relative distance between nodes in east direction
nodes_north	relative distance between nodes in north direction
nodes_z	relative distance between nodes in east direction
pad_east	number of cells for padding on E and W sides default is 7
pad_north	number of cells for padding on S and N sides default is 7
pad_num	number of cells with cell_size with outside of station area. default is 3
pad_method	method to use to create padding: extent1, extent2 - calculate based on ew_ext and ns_ext stretch - calculate based on pad_stretch factors
pad_stretch_h	multiplicative number for padding in horizontal direction.
pad_stretch_v	padding cells N & S will be pad_root_north**(x)
pad_z	number of cells for padding at bottom default is 4
ew_ext	E-W extension of model in meters
ns_ext	N-S extension of model in meters
res_scale	scaling method of res, supports ‘loge’ - for log e format ‘log’ or ‘log10’ - for log with base 10 ‘linear’ - linear scale default is ‘loge’
res_list	list of resistivity values for starting model
res_model	starting resistivity model
res_initial_value	resistivity initial value for the resistivity model default is 100
mesh_rotation_angle	Angle to rotate the grid to. Angle is measured positve clockwise assuming North is 0 and east is 90. default is None
save_path	path to save file to
sea_level	sea level in grid_z coordinates. default is 0
station_locations	location of stations
title	title in initial file
z1_layer	first layer thickness
z_bottom	absolute bottom of the model default is 300,000
z_target_depth	Depth of deepest target, default is 50,000

add_layers_to_mesh(n_add_layers=None, layer_thickness=None, where='top')[source]

Function to add constant thickness layers to the top or bottom of mesh.

Note: It is assumed these layers are added before the topography. If you want to add topography layers, use function add_topography_to_model :param n_add_layers: Integer, number of layers to add, defaults to None. :param layer_thickness: Real value or list/array. Thickness of layers,

Can provide a single value

or a list/array containing multiple layer thicknesses, defaults to None.

Parameters:: where – Where to add, top or bottom, defaults to “top”.

add_topography_from_data(interp_method='nearest', air_resistivity=1000000000000.0, topography_buffer=None, airlayer_type='log_up')[source]

Wrapper around add_topography_to_model that allows creating a surface model from EDI data. The Data grid and station elevations will be used to make a ‘surface’ tuple that will be passed to add_topography_to_model so a surface model can be interpolated from it.

The surface tuple is of format (lon, lat, elev) containing station locations. :param data_object: A ModEm data

object that has been filled with data from EDI files.

Parameters:

interp_method (str, optional) – Same as add_topography_to_model, defaults to “nearest”.
air_resistivity (float, optional) – Same as add_topography_to_model, defaults to 1e12.
topography_buffer (float, optional) – Same as add_topography_to_model, defaults to None.
airlayer_type (str, optional) – Same as add_topography_to_model, defaults to “log_up”.

add_topography_to_model(topography_file=None, surface=None, topography_array=None, interp_method='nearest', air_resistivity=1000000000000.0, topography_buffer=None, airlayer_type='log_up', max_elev=None, shift_east=0, shift_north=0)[source]

If air_layers is non-zero, will add topo: read in topograph file, make a surface model.

Call project_stations_on_topography in the end, which will re-write the .dat file.

If n_airlayers is zero, then cannot add topo data, only bathymetry is needed. :param shift_north:

Defaults to 0.

Parameters:

shift_east – Defaults to 0.
max_elev – Defaults to None.
surface – Defaults to None.
topography_file – File containing topography (arcgis ascii grid), defaults to None.
topography_array – Alternative to topography_file - array of elevation values on model grid, defaults to None.
interp_method – Interpolation method for topography, ‘nearest’, ‘linear’, or ‘cubic’, defaults to “nearest”.
air_resistivity – Resistivity value to assign to air, defaults to 1e12.
topography_buffer – Buffer around stations to calculate minimum and maximum topography value to use for meshing, defaults to None.
airlayer_type – How to set air layer thickness - options are ‘constant’ for constant air layer thickness, or ‘log’, for logarithmically increasing air layer thickness upward, defaults to “log_up”.

assign_resistivity_from_surface_data(top_surface, bottom_surface, resistivity_value)[source]

Assign resistivity value to all points above or below a surface requires the surface_dict attribute to exist and contain data for surface key (can get this information from ascii file using project_surface)

inputs surface_name = name of surface (must correspond to key in surface_dict) resistivity_value = value to assign where = ‘above’ or ‘below’ - assign resistivity above or below the

surface

interpolate_elevation(surface_file=None, surface=None, get_surface_name=False, method='nearest', fast=True, shift_north=0, shift_east=0)[source]

Project a surface to the model grid and add resulting elevation data to a dictionary called surface_dict. Assumes the surface is in lat/long coordinates (wgs84)

returns nothing returned, but surface data are added to surface_dict under the key given by surface_name.

inputs choose to provide either surface_file (path to file) or surface (tuple). If both are provided then surface tuple takes priority.

surface elevations are positive up, and relative to sea level. surface file format is:

ncols 3601 nrows 3601 xllcorner -119.00013888889 (longitude of lower left) yllcorner 36.999861111111 (latitude of lower left) cellsize 0.00027777777777778 NODATA_value -9999 elevation data W –> E N | V S

Alternatively, provide a tuple with: (lon,lat,elevation) where elevation is a 2D array (shape (ny,nx)) containing elevation points (order S -> N, W -> E) and lon, lat are either 1D arrays containing list of longitudes and latitudes (in the case of a regular grid) or 2D arrays with same shape as elevation array containing longitude and latitude of each point.

other inputs: surface_epsg = epsg number of input surface, default is 4326 for lat/lon(wgs84) method = interpolation method. Default is ‘nearest’, if model grid is dense compared to surface points then choose ‘linear’ or ‘cubic’

make_mesh(verbose=True)[source]

Create finite element mesh according to user-input parameters.

The mesh is built by:

Making a regular grid within the station area.
Adding pad_num of cell_width cells outside of station area
Adding padding cells to given extension and number of padding cells.
Making vertical cells starting with z1_layer increasing logarithmically (base 10) to z_target_depth and num_layers.
Add vertical padding cells to desired extension.
Check to make sure none of the stations lie on a node. If they do then move the node by .02*cell_width

make_z_mesh(n_layers=None)[source]: New version of make_z_mesh. make_z_mesh and M.

property model_epsg: Model epsg.

property model_fn: Model fn.

property model_parameters: Get important model parameters to write to a file for documentation later.

property nodes_east: Nodes east.

property nodes_north: Nodes north.

property nodes_z: Nodes z.

property plot_east: Plot east.

plot_mesh(**kwargs)[source]: Plot model mesh. :param **kwargs: :param plot_topography: DESCRIPTION, defaults to False. :type plot_topography: TYPE, optional :return: DESCRIPTION. :rtype: TYPE

property plot_north: Plot north.

property plot_z: Plot z.

read_gocad_sgrid_file(sgrid_header_file, air_resistivity=1e+39, sea_resistivity=0.3, sgrid_positive_up=True)[source]

Read a gocad sgrid file and put this info into a ModEM file.

Note: can only deal with grids oriented N-S or E-W at this stage, with orthogonal coordinates

read_model_file(model_fn=None)[source]

Read an initial file and return the pertinent information including grid positions in coordinates relative to the center point (0,0) and starting model.

Note that the way the model file is output, it seems is that the blocks are setup as

ModEM: WS::: 0—–> N_north 0——–>N_east | | | | V V N_east N_north

Arguments:

**model_fn** : full path to initializing file.

Outputs:

**nodes_north** : np.array(nx)
            array of nodes in S --> N direction

**nodes_east** : np.array(ny)
            array of nodes in the W --> E direction

**nodes_z** : np.array(nz)
            array of nodes in vertical direction positive downwards

**res_model** : dictionary
            dictionary of the starting model with keys as layers

**res_list** : list
            list of resistivity values in the model

**title** : string
             title string

property save_path: Save path.

write_geosoft_xyz(save_fn, c_east=0, c_north=0, c_z=0, pad_north=0, pad_east=0, pad_z=0)[source]

Write an XYZ file readable by Geosoft

All input units are in meters.. :param save_fn: Full path to save file to. :type save_fn: string or Path :param c_east: Center point in the east direction, defaults to 0. :type c_east: float, optional :param c_north: Center point in the north direction, defaults to 0. :type c_north: float, optional :param c_z: Center point elevation, defaults to 0. :type c_z: float, optional :param pad_north: Number of cells to cut from the north-south edges, defaults to 0. :type pad_north: int, optional :param pad_east: Number of cells to cut from the east-west edges, defaults to 0. :type pad_east: int, optional :param pad_z: Number of cells to cut from the bottom, defaults to 0. :type pad_z: int, optional

write_gocad_sgrid_file(fn=None, origin=[0, 0, 0], clip=0, no_data_value=-99999)[source]

Write a model to gocad sgrid

optional inputs:

fn = filename to save to. File extension (‘.sg’) will be appended.: default is the model name with extension removed

origin = real world [x,y,z] location of zero point in model grid clip = how much padding to clip off the edge of the model for export,

provide one integer value or list of 3 integers for x,y,z directions

no_data_value = no data value to put in sgrid.

write_model_file(**kwargs)[source]

Will write an initial file for ModEM.

Note that x is assumed to be S –> N, y is assumed to be W –> E and z is positive downwards. This means that index [0, 0, 0] is the southwest corner of the first layer. Therefore if you build a model by hand the layer block will look as it should in map view.

Also, the xgrid, ygrid and zgrid are assumed to be the relative distance between neighboring nodes. This is needed because wsinv3d builds the model from the bottom SW corner assuming the cell width from the init file.

Key Word Arguments:

**nodes_north** : np.array(nx)
            block dimensions (m) in the N-S direction.
            **Note** that the code reads the grid assuming that
            index=0 is the southern most point.

**nodes_east** : np.array(ny)
            block dimensions (m) in the E-W direction.
            **Note** that the code reads in the grid assuming that
            index=0 is the western most point.

**nodes_z** : np.array(nz)
            block dimensions (m) in the vertical direction.
            This is positive downwards.

**save_path** : string
              Path to where the initial file will be saved
              to save_path/model_fn_basename

**model_fn_basename** : string
                        basename to save file to
                        *default* is ModEM_Model.ws
                        file is saved at save_path/model_fn_basename

**title** : string
            Title that goes into the first line
            *default* is Model File written by MTpy.modeling.modem

**res_model** : np.array((nx,ny,nz))
            Prior resistivity model.

            .. note:: again that the modeling code
            assumes that the first row it reads in is the southern
            most row and the first column it reads in is the
            western most column.  Similarly, the first plane it
            reads in is the Earth's surface.

**res_starting_value** : float
                         starting model resistivity value,
                         assumes a half space in Ohm-m
                         *default* is 100 Ohm-m

**res_scale** : [ 'loge' | 'log' | 'log10' | 'linear' ]
                scale of resistivity.  In the ModEM code it
                converts everything to Loge,
                *default* is 'loge'

write_out_file(save_fn, geographic_east, geographic_north, geographic_elevation)[source]

Will write an .out file for LeapFrog.

Note that y is assumed to be S –> N, e is assumed to be W –> E and z is positive upwards. This means that index [0, 0, 0] is the southwest corner of the first layer. :param save_fn: Full path to save file to. :type save_fn: string or Path :param geographic_east: Geographic center in easting (meters). :type geographic_east: float :param geographic_north: Geographic center in northing (meters). :type geographic_north: float :param geographic_elevation: Elevation of geographic center (meters). :type geographic_elevation: float :return: DESCRIPTION. :rtype: TYPE

write_ubc_files(basename, c_east=0, c_north=0, c_z=0)[source]: Write a UBC .msh and .mod file. :param basename: :param save_fn: DESCRIPTION. :type save_fn: TYPE :param c_east: DESCRIPTION, defaults to 0. :type c_east: TYPE, optional :param c_north: DESCRIPTION, defaults to 0. :type c_north: TYPE, optional :param c_z: DESCRIPTION, defaults to 0. :type c_z: TYPE, optional :return: DESCRIPTION. :rtype: TYPE

write_vtk_file(vtk_save_path=None, vtk_fn_basename='ModEM_model_res', shift_east=0, shift_north=0, shift_z=0, units='km', coordinate_system='nez+', label='resistivity')[source]

Write a VTK file to plot in 3D rendering programs like Paraview. :param label:

Defaults to “resistivity”.

Parameters:

coordinate_system – Defaults to “nez+”.
units – Defaults to “km”.
shift_z – Defaults to 0.
shift_north – Defaults to 0.
shift_east – Defaults to 0.
vtk_save_path (string or Path, optional) – Directory to save vtk file to, defaults to None.
vtk_fn_basename – Filename basename of vtk file, note that .vtr, defaults to “ModEM_model_res”.

Returns:

Full path to VTK file.

Return type:

Path

write_xyres(save_path=None, location_type='EN', origin=[0, 0], model_epsg=None, depth_index='all', outfile_basename='DepthSlice', log_res=False, clip=[0, 0])[source]

Write files containing depth slice data (x, y, res for each depth)

origin = x,y coordinate of zero point of ModEM_grid, or name of file: containing this info (full path or relative to model files)

save_path = path to save to, default is the model object save path location_type = ‘EN’ or ‘LL’ xy points saved as eastings/northings or

longitude/latitude, if ‘LL’ need to also provide model_epsg

model_epsg = epsg number that was used to project the model outfile_basename = string for basename for saving the depth slices. log_res = True/False - option to save resistivity values as log10

instead of linear

clip = number of cells to clip on each of the east/west and north/south edges.

write_xyzres(savefile=None, location_type='EN', origin=[0, 0], model_epsg=None, log_res=False, model_utm_zone=None, clip=[0, 0])[source]: Save a model file as a space delimited x y z res file.

mtpy.modeling.modem.residual module

ModEM

residuals class to contain RMS information

revised by JP 2017 revised by AK 2017 to bring across functionality from ak branch

class mtpy.modeling.modem.residual.Residual(**kwargs)[source]

Bases: Data

Class to contain residuals for each data point, and rms values for each station

Attributes/Key Words	Description
work_dir
residual_fn	full path to data file
residual_array	numpy.ndarray (num_stations) structured to store data. keys are: station –> station name lat –> latitude in decimal degrees lon –> longitude in decimal degrees elev –> elevation (m) rel_east – > relative east location to center_position (m) rel_north –> relative north location to center_position (m) east –> UTM east (m) north –> UTM north (m) zone –> UTM zone z –> impedance tensor residual (measured - modelled) (num_freq, 2, 2) z_err –> impedance tensor error array with shape (num_freq, 2, 2) tip –> Tipper residual (measured - modelled) (num_freq, 1, 2) tipperr –> Tipper array with shape (num_freq, 1, 2)
rms
rms_array	numpy.ndarray structured to store station location values and rms. Keys are: station –> station name east –> UTM east (m) north –> UTM north (m) lat –> latitude in decimal degrees lon –> longitude in decimal degrees elev –> elevation (m) zone –> UTM zone rel_east – > relative east location to center_position (m) rel_north –> relative north location to center_position (m) rms –> root-mean-square residual for each station
rms_tip
rms_z

calculate_rms()[source]: Add columns for rms. :return: DESCRIPTION. :rtype: TYPE

plot_rms(**kwargs)[source]: Plot RMS in different views. :param **kwargs: DESCRIPTION. :type **kwargs: TYPE :return: DESCRIPTION. :rtype: TYPE

plot_rms_per_period(plot_type='all', **kwargs)[source]

Plot rms per period. :param plot_type:

Defaults to “all”.

Parameters:

**kwargs –

DESCRIPTION.

Returns:

DESCRIPTION.

Return type:

TYPE

read_residual_file(residual_fn)[source]: Read residual file. :param residual_fn: DESCRIPTION, defaults to None. :type residual_fn: TYPE, optional :return: DESCRIPTION. :rtype: TYPE

property rms_per_period_all: RMS per period.

property rms_per_period_per_component: RMS per period by component. :return: DESCRIPTION. :rtype: TYPE

mtpy.modeling.modem.station module

ModEM

# Generate files for ModEM

# revised by JP 2017 # revised by AK 2017 to bring across functionality from ak branch

class mtpy.modeling.modem.station.Stations(**kwargs)[source]

Bases: object

Station locations class.

calculate_rel_locations(shift_east=0, shift_north=0)[source]: Put station in a coordinate system relative to (shift_east, shift_north) (+) shift right or up (-) shift left or down

property center_point: Calculate the center point from the given station locations.

property east: East function.

property elev: Elev function.

get_station_locations(input_list)[source]

Get station locations from a list of edi files.

Parameters:: **input_list** – list list of edi file names, or mt_objects

property lat: Lat function.

property lon: Lon function.

property model_epsg: Model epsg.

property model_utm_zone: Model utm zone.

property north: North function.

property rel_east: Rel east.

property rel_elev: Rel elev.

property rel_north: Rel north.

rotate_stations(rotation_angle)[source]

Rotate stations assuming N is 0.

Parameters:: **rotation_angle** – float angle in degrees assuming N is 0

property station: Station function.

to_csv(csv_fn, epsg=None, default_epsg=4326, geometry=False)[source]

Write a shape file of the station locations using geopandas which only takes in epsg numbers :param geometry:

Defaults to False.

Parameters:

csv_fn
shp_fn (string) – Full path to new shapefile.
epsg (integer, defaults to None, optional) – EPSG number to project to, defaults to None.
default_epsg – The default EPSG number that the stations are, defaults to 4326.

to_geopd(epsg=None, default_epsg=4326)[source]: Create a geopandas dataframe. :param epsg: EPSG number to project to, defaults to None. :type epsg: integer, defaults to None, optional :param default_epsg: The default EPSG number that the stations are, defaults to 4326.

to_shp(shp_fn, epsg=None, default_epsg=4326)[source]: Write a shape file of the station locations using geopandas which only takes in epsg numbers :param shp_fn: Full path to new shapefile. :type shp_fn: string :param epsg: EPSG number to project to, defaults to None. :type epsg: integer, defaults to None, optional :param default_epsg: The default EPSG number that the stations are, defaults to 4326.

property utm_zone: Utm zone.

Module contents

class mtpy.modeling.modem.ControlFwd(**kwargs)[source]

Bases: object

Reads and writes control files for ModEM to control how the forward solver starts and how it is run.

num_qmr_iter Number of QMR iterations per divergence correction.: int, default is 40
max_num_div_calls Maximum number of divergence correction calls. int,: default is 20
max_num_div_iters Maximum number of divergence correction iterations.: int, default is 100
misfit_tol_fwd Misfit tolerance for EM forward solver. float,: default is 1e-07
misfit_tol_adj Misfit tolerance for EM adjoint solver. float,: default is 1e-07
misfit_tol_div Misfit tolerance for divergence correction. float,: default is 1e-05
save_path Path at which to save the control file. PosixPath,: default is current working directory.
fn_basename Name of the control file. str, default is: ‘control.fwd’,

property control_fn: Control fn.

read_control_file(control_fn=None)[source]: Read in a control file.

write_control_file(control_fn=None, save_path=None, fn_basename=None)[source]

Write control file.

Arguments::

control_fnstring: full path to save control file to default is save_path/fn_basename
save_pathstring: directory path to save control file to default is cwd
fn_basenamestring: basename of control file default is control.inv

class mtpy.modeling.modem.ControlInv(**kwargs)[source]

Bases: object

Reads and writes control files for ModEM to control how the inversion starts and how it is run.

output_fn Model and data output file name “prefix phrase”, i.e.: written to the front of ModEM output file names before the iteration number and file extension. str, default is ‘MODULAR_NLCG’

lambda_initial Initial damping factor lambda. int, default is 10 lambda_step To update lambda, divide by this value. int, default

is 10

model_search_step Initial search step in model units. int, default is: 1
rms_reset_search Restart when rms diff is less than this value. float,: default is 0.002
rms_target Exit search when rms is less than this value. float,: default is 1.05
lambda_exit Exit when lambda is less than this value. float,: default is 0.0001

max_iterations Maximum number of iterations. int, default is 100 save_path Path at which to save the control file. PosixPath,

default is current working directory

fn_basename Name of the control file. str, default is: ‘control.inv’

property control_fn: Control fn.

read_control_file(control_fn=None)[source]: Read in a control file.

write_control_file(control_fn=None, save_path=None, fn_basename=None)[source]

Write control file.

Arguments::

control_fnstring: full path to save control file to default is save_path/fn_basename
save_pathstring: directory path to save control file to default is cwd
fn_basenamestring: basename of control file default is control.inv

class mtpy.modeling.modem.Covariance(grid_dimensions=None, **kwargs)[source]

Bases: object

Class that deals with model covariance and smoothing, writing covariance (.cov) files for use in ModEM.

grid_dimensions Grid dimensions excluding air layers (Nx, Ny, NzEarth) smoothing_east Smoothing in east direction. float, default is 0.3 smoothing_north Smoothing in north direction. float, default is 0.3 smoothing_z Smoothing in vertical direction. float, default is

0.3

smoothing_num Number of smoothing iterations. int, default is 1 exception_list List of exceptions. list, default is empty mask_arr Mask array. numpy.ndarray, default is None save_path Path at which to save the covariance file. PosixPath,

default is current working directory

fn_basename Name of the covariance file. default is: ‘covariance.cov’

property cov_fn: Cov fn.

get_parameters()[source]: Get parameters.

read_cov_file(cov_fn)[source]: Reads a ModEM covariance (.cov) file. :param cov_fn: Filename of the target ModEM covariance (.cov) file. :type cov_fn: str :raises CovarianceError: Error related to the covariance class.

write_cov_vtk_file(cov_vtk_fn, model_fn=None, grid_east=None, grid_north=None, grid_z=None)[source]: Write a vtk file of the covariance to match things up.

write_covariance_file(cov_fn=None, save_path=None, fn_basename=None, res_model=None, sea_water=0.3, air=1000000000000.0)[source]

Writes a ModEM covariance (.cov) file. :param cov_fn: Name for the covariance file. default is None, defaults to None. :type cov_fn: str | Path, optional :param save_path: Location at which to save the covariance file.

default is None, defaults to None.

Parameters:

fn_basename (_type_, optional) – _description_, default is None, defaults to None.
res_model (_type_, optional) – Resistivity model the covariance should be generated from, default is None, defaults to None.
sea_water (float, optional) – Electrical resistivity of sea water cells in Ohm-meters within the model domain, default is 0.3, defaults to 0.3.
air (float, optional) – Electrical resistivity of air cells in Ohm-meters within the model domain, default is 1.0e12, defaults to 1.0e12.

Raises:

CovarianceError – Error related to the covariance class.

class mtpy.modeling.modem.Data(dataframe=None, center_point=None, **kwargs)[source]

Bases: object

Data will read and write .dat files for ModEM and convert a WS data file to ModEM format.

..note: :: the data is interpolated onto the given periods such that all: stations invert for the same periods. The interpolation is a linear interpolation of each of the real and imaginary parts of the impedance tensor and induction tensor. See mtpy.core.mt.MT.interpolate for more details

Parameters:

**kwargs –
center_point – Defaults to None.
dataframe – Defaults to None.
edi_list – List of edi files to read.

property data_filename: Data filename.

property dataframe: Dataframe function.

fix_data_file(fn=None, n=3)[source]: A newer compiled version of Modem outputs slightly different headers This aims to convert that into the older format :param fn: DESCRIPTION, defaults to None. :type fn: TYPE, optional :param n: DESCRIPTION, defaults to 3. :type n: TYPE, optional :return: DESCRIPTION. :rtype: TYPE

get_n_stations(mode)[source]: Get n stations.

property model_parameters: Model parameters.

property n_periods: N periods.

property period: Period function.

read_data_file(data_fn)[source]: Read data file. :param data_fn: Full path to data file name. :type data_fn: string or Path :raises ValueError: If cannot compute component.

write_data_file(file_name=None, save_path=None, fn_basename=None, elevation=False)[source]

Write data file. :param file_name:

Defaults to None.

Parameters:

save_path (string or Path, optional) – Full directory to save file to, defaults to None.
fn_basename (string, optional) – Basename of the saved file, defaults to None.
elevation (boolean, optional) – If True adds in elevation from ‘rel_elev’ column in data array, defaults to False.

Raises:

NotImplementedError – If the inversion mode is not supported.
ValueError – mtpy.utils.exceptions.ValueError if a parameter.

Returns:

Full path to data file.

Return type:

Path

exception mtpy.modeling.modem.DataError[source]

Bases: ModEMError

Raise for ModEM Data class specific exceptions.

class mtpy.modeling.modem.ModEMConfig(**kwargs)[source]

Bases: object

Read and write configuration files for how each inversion is run.

add_dict(fn=None, obj=None)[source]: Add dictionary based on file name or object.

write_config_file(save_dir=None, config_fn_basename='ModEM_inv.cfg')[source]: Write a config file based on provided information.

exception mtpy.modeling.modem.ModEMError[source]: Bases: Exception

class mtpy.modeling.modem.Model(station_locations=None, center_point=None, **kwargs)[source]

Bases: object

Make and read a FE mesh grid

The mesh assumes the coordinate system where:: x == North y == East z == + down

All dimensions are in meters.

The mesh is created by first making a regular grid around the station area, then padding cells are added that exponentially increase to the given extensions. Depth cell increase on a log10 scale to the desired depth, then padding cells are added that increase exponentially..

Parameters:: **station_object** – mtpy.modeling.modem.Stations object .. seealso:: mtpy.modeling.modem.Stations

Examples

Example 1 –> create mesh first then data file:

>>> import mtpy.modeling.modem as modem
>>> import os
>>> # 1) make a list of all .edi files that will be inverted for
>>> edi_path = r"/home/EDI_Files"
>>> edi_list = [os.path.join(edi_path, edi)

for edi in os.listdir(edi_path)

>>> ...         if edi.find('.edi') > 0]
>>> # 2) Make a Stations object
>>> stations_obj = modem.Stations()
>>> stations_obj.get_station_locations_from_edi(edi_list)
>>> # 3) make a grid from the stations themselves with 200m cell spacing
>>> mmesh = modem.Model(station_obj)
>>> # change cell sizes
>>> mmesh.cell_size_east = 200,
>>> mmesh.cell_size_north = 200
>>> mmesh.ns_ext = 300000 # north-south extension
>>> mmesh.ew_ext = 200000 # east-west extension of model
>>> mmesh.make_mesh()
>>> # check to see if the mesh is what you think it should be
>>> msmesh.plot_mesh()
>>> # all is good write the mesh file
>>> msmesh.write_model_file(save_path=r"/home/modem/Inv1")
>>> # create data file
>>> md = modem.Data(edi_list, station_locations=mmesh.station_locations)
>>> md.write_data_file(save_path=r"/home/modem/Inv1")

Example 2 –> Rotate Mesh:

>>> mmesh.mesh_rotation_angle = 60
>>> mmesh.make_mesh()

Note

ModEM assumes all coordinates are relative to North and East, and does not accommodate mesh rotations, therefore, here the rotation is of the stations, which essentially does the same thing. You will need to rotate you data to align with the ‘new’ coordinate system.

Attributes	Description
_logger	python logging object that put messages in logging format defined in logging configure file, see MtPyLog more information
cell_number_ew	optional for user to specify the total number of sells on the east-west direction. default is None
cell_number_ns	optional for user to specify the total number of sells on the north-south direction. default is None
cell_size_east	mesh block width in east direction default is 500
cell_size_north	mesh block width in north direction default is 500
grid_center	center of the mesh grid
grid_east	overall distance of grid nodes in east direction
grid_north	overall distance of grid nodes in north direction
grid_z	overall distance of grid nodes in z direction
model_fn	full path to initial file name
model_fn_basename	default name for the model file name
n_air_layers	number of air layers in the model. default is 0
n_layers	total number of vertical layers in model
nodes_east	relative distance between nodes in east direction
nodes_north	relative distance between nodes in north direction
nodes_z	relative distance between nodes in east direction
pad_east	number of cells for padding on E and W sides default is 7
pad_north	number of cells for padding on S and N sides default is 7
pad_num	number of cells with cell_size with outside of station area. default is 3
pad_method	method to use to create padding: extent1, extent2 - calculate based on ew_ext and ns_ext stretch - calculate based on pad_stretch factors
pad_stretch_h	multiplicative number for padding in horizontal direction.
pad_stretch_v	padding cells N & S will be pad_root_north**(x)
pad_z	number of cells for padding at bottom default is 4
ew_ext	E-W extension of model in meters
ns_ext	N-S extension of model in meters
res_scale	scaling method of res, supports ‘loge’ - for log e format ‘log’ or ‘log10’ - for log with base 10 ‘linear’ - linear scale default is ‘loge’
res_list	list of resistivity values for starting model
res_model	starting resistivity model
res_initial_value	resistivity initial value for the resistivity model default is 100
mesh_rotation_angle	Angle to rotate the grid to. Angle is measured positve clockwise assuming North is 0 and east is 90. default is None
save_path	path to save file to
sea_level	sea level in grid_z coordinates. default is 0
station_locations	location of stations
title	title in initial file
z1_layer	first layer thickness
z_bottom	absolute bottom of the model default is 300,000
z_target_depth	Depth of deepest target, default is 50,000

add_layers_to_mesh(n_add_layers=None, layer_thickness=None, where='top')[source]

Function to add constant thickness layers to the top or bottom of mesh.

Note: It is assumed these layers are added before the topography. If you want to add topography layers, use function add_topography_to_model :param n_add_layers: Integer, number of layers to add, defaults to None. :param layer_thickness: Real value or list/array. Thickness of layers,

Can provide a single value

or a list/array containing multiple layer thicknesses, defaults to None.

Parameters:: where – Where to add, top or bottom, defaults to “top”.

add_topography_from_data(interp_method='nearest', air_resistivity=1000000000000.0, topography_buffer=None, airlayer_type='log_up')[source]

Wrapper around add_topography_to_model that allows creating a surface model from EDI data. The Data grid and station elevations will be used to make a ‘surface’ tuple that will be passed to add_topography_to_model so a surface model can be interpolated from it.

The surface tuple is of format (lon, lat, elev) containing station locations. :param data_object: A ModEm data

object that has been filled with data from EDI files.

Parameters:

interp_method (str, optional) – Same as add_topography_to_model, defaults to “nearest”.
air_resistivity (float, optional) – Same as add_topography_to_model, defaults to 1e12.
topography_buffer (float, optional) – Same as add_topography_to_model, defaults to None.
airlayer_type (str, optional) – Same as add_topography_to_model, defaults to “log_up”.

add_topography_to_model(topography_file=None, surface=None, topography_array=None, interp_method='nearest', air_resistivity=1000000000000.0, topography_buffer=None, airlayer_type='log_up', max_elev=None, shift_east=0, shift_north=0)[source]

If air_layers is non-zero, will add topo: read in topograph file, make a surface model.

Call project_stations_on_topography in the end, which will re-write the .dat file.

If n_airlayers is zero, then cannot add topo data, only bathymetry is needed. :param shift_north:

Defaults to 0.

Parameters:

shift_east – Defaults to 0.
max_elev – Defaults to None.
surface – Defaults to None.
topography_file – File containing topography (arcgis ascii grid), defaults to None.
topography_array – Alternative to topography_file - array of elevation values on model grid, defaults to None.
interp_method – Interpolation method for topography, ‘nearest’, ‘linear’, or ‘cubic’, defaults to “nearest”.
air_resistivity – Resistivity value to assign to air, defaults to 1e12.
topography_buffer – Buffer around stations to calculate minimum and maximum topography value to use for meshing, defaults to None.
airlayer_type – How to set air layer thickness - options are ‘constant’ for constant air layer thickness, or ‘log’, for logarithmically increasing air layer thickness upward, defaults to “log_up”.