mtpy.modeling package

Subpackages

Submodules

mtpy.modeling.errors module

Created on Tue Oct 11 16:01:37 2022

@author: jpeacock

class mtpy.modeling.errors.ModelErrors(data=None, measurement_error=None, **kwargs)[source]

Bases: object

compute_absolute_error()[source]: Compute absolute error. :param data: DESCRIPTION. :type data: TYPE :param error_value: DESCRIPTION. :type error_value: TYPE :return: DESCRIPTION. :rtype: TYPE

compute_arithmetic_mean_error()[source]: Error_value * (Zxy + Zyx) / 2. :param data: DESCRIPTION. :type data: TYPE :param error_value: DESCRIPTION. :type error_value: TYPE :return: DESCRIPTION. :rtype: TYPE

compute_eigen_value_error()[source]: Error_value * eigen(data).mean(). :param data: DESCRIPTION. :type data: TYPE :param error_value: DESCRIPTION. :type error_value: TYPE :return: DESCRIPTION. :rtype: TYPE

compute_error(data=None, error_type=None, error_value=None, floor=None)[source]: Compute error. :param data: DESCRIPTION, defaults to None. :type data: TYPE, optional :param error_type: DESCRIPTION, defaults to None. :type error_type: TYPE, optional :param error_value: DESCRIPTION, defaults to None. :type error_value: TYPE, optional :param floor: DESCRIPTION, defaults to None. :type floor: TYPE, optional :return: DESCRIPTION. :rtype: TYPE

compute_geometric_mean_error()[source]: Error_value * sqrt(Zxy * Zyx). :param data: DESCRIPTION. :type data: TYPE :param error_value: DESCRIPTION. :type error_value: TYPE :return: DESCRIPTION. :rtype: TYPE

compute_median_error()[source]: Median(array) * error_value. :param array: DESCRIPTION. :type array: TYPE :param error_value: DESCRIPTION. :type error_value: TYPE :return: DESCRIPTION. :rtype: TYPE

compute_percent_error()[source]: Percent error. :param data: DESCRIPTION. :type data: TYPE :param percent: DESCRIPTION. :type percent: TYPE :return: DESCRIPTION. :rtype: TYPE

compute_row_error()[source]: Set zxx and zxy the same error and zyy and zyx the same error. :param data: DESCRIPTION. :type data: TYPE :param error_value: DESCRIPTION. :type error_value: TYPE :param floor: DESCRIPTION, defaults to True. :type floor: TYPE, optional :return: DESCRIPTION. :rtype: TYPE

property data: Data function.

property error_parameters: Error parameters.

property error_type: Error type.

property error_value: Error value.

property floor: Floor function.

mask_zeros(data)[source]: Mask zeros. :param data: DESCRIPTION. :type data: TYPE :return: DESCRIPTION. :rtype: TYPE

property measurement_error: Measurement error.

property mode: Mode function.

resize_output(error_array)[source]: Resize the error estimtion to the same size as the input data. :param error_array: DESCRIPTION. :type error_array: TYPE :return: DESCRIPTION. :rtype: TYPE

set_floor(error_array)[source]: Set error floor. :param error_array: :param array: DESCRIPTION. :type array: TYPE :param floor: DESCRIPTION. :type floor: TYPE :return: DESCRIPTION. :rtype: TYPE

use_measurement_error()[source]: Use measurement error.

validate_array_shape(data)[source]: Validate array shape. :param data: DESCRIPTION. :type data: TYPE :return: DESCRIPTION. :rtype: TYPE

validate_percent(value)[source]: Make sure the percent is a decimal. :param value: DESCRIPTION. :type value: TYPE :return: DESCRIPTION. :rtype: TYPE

mtpy.modeling.gocad module

Created on Fri Dec 09 15:50:53 2016

@author: Alison Kirkby

read and write gocad objects

class mtpy.modeling.gocad.Sgrid(**kwargs)[source]

Bases: object

Class to read and write gocad sgrid files

need to provide: workdir = working directory fn = filename for the sgrid resistivity = 3d numpy array containing resistivity values, shape (ny,nx,nz) grid_xyz = tuple containing x,y,z locations of edges of cells for each

resistivity value. Each item in tuple has shape (ny+1,nx+1,nz+1).

read_sgrid_file(headerfn=None)[source]: Read sgrid file.

write_sgrid_file()[source]: Write sgrid file.

mtpy.modeling.mare2dem module

mtpy.modeling.mesh_tools module

Created on Wed Oct 25 09:35:31 2017

@author: Alison Kirkby

functions to assist with mesh generation

mtpy.modeling.mesh_tools.get_nearest_index(array, value)[source]

Return the index of the nearest value to the provided value in an array:

inputs:: array = array or list of values value = target value.

mtpy.modeling.mesh_tools.get_padding_cells(cell_width, max_distance, num_cells, stretch)[source]

Get padding cells, which are exponentially increasing to a given distance. Make sure that each cell is larger than the one previously.

Parameters:

**cell_width** – float width of grid cell (m)
**max_distance** – float maximum distance the grid will extend (m)
**num_cells** – int number of padding cells
**stretch** – float base geometric factor

mtpy.modeling.mesh_tools.get_padding_cells2(cell_width, core_max, max_distance, num_cells)[source]: Get padding cells, which are exponentially increasing to a given distance. Make sure that each cell is larger than the one previously.

mtpy.modeling.mesh_tools.get_padding_from_stretch(cell_width, pad_stretch, num_cells)[source]: Get padding cells using pad stretch factor.

mtpy.modeling.mesh_tools.get_rounding(cell_width)[source]

Get the rounding number given the cell width.

Will be one significant number less

than the cell width. This reduces weird looking meshes.

param cell_width:: Width of mesh cell.
type cell_width:: float
return:: Digit to round to.
rtype:: int

mtpy.modeling.mesh_tools.get_station_buffer(grid_east, grid_north, station_east, station_north, buf=10000.0)[source]: Get cells within a specified distance (buf) of the stations returns a 2D boolean (True/False) array

mtpy.modeling.mesh_tools.grid_centre(grid_edges)[source]: Calculate the grid centres from an array that defines grid edges. :param grid_edges: Array containing grid edges. :return s: Grid_centre: centre points of grid.

mtpy.modeling.mesh_tools.interpolate_elevation_to_grid(grid_east, grid_north, utm_epsg=None, datum_epsg=4326, surface_file=None, surface=None, method='linear', fast=True, buffer=1)[source]

# Note: this documentation is outdated and seems to be copied from # model.interpolate_elevation2. It needs to be updated. This # funciton does not update a dictionary but returns an array of # elevation data.

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) The ‘fast’ method extracts a subset of the elevation data that falls within the mesh-bounds and interpolates them onto mesh nodes. This approach significantly speeds up (~ x5) the interpolation procedure.

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

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: surfacename = name of surface for putting into dictionary 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’

mtpy.modeling.mesh_tools.make_log_increasing_array(z1_layer, target_depth, n_layers, increment_factor=0.9)[source]: Create depth array with log increasing cells, down to target depth, inputs are z1_layer thickness, target depth, number of layers (n_layers)

mtpy.modeling.mesh_tools.rotate_mesh(grid_east, grid_north, origin, rotation_angle, return_centre=False)[source]: Rotate a mesh defined by grid_east and grid_north. :param grid_east: 1d array defining the edges of the mesh in the east-west direction. :param grid_north: 1d array defining the edges of the mesh in the north-south direction. :param origin: Real-world position of the (0,0) point in grid_east, grid_north. :param rotation_angle: Angle in degrees to rotate the grid by. :param return_centre: True/False option to return points on centre of grid instead of grid edges, defaults to False. :return: Grid_east, grid_north - 2d arrays describing the east and north coordinates.

mtpy.modeling.pek1d module

mtpy.modeling.pek1dclasses module

mtpy.modeling.pek2d module

mtpy.modeling.pek2dforward module

mtpy.modeling.structured_mesh_3d 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.structured_mesh_3d.StructuredGrid3D(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

convert_model_to_int(res_list=None)[source]: Convert resistivity values to integers according to resistivity list. :param res_list: Resistivity values in Ohm-m, defaults to None. :type res_list: list of floats, optional :return: Array of integers corresponding to the res_list. :rtype: np.ndarray(dtype=int)

estimate_skin_depth(apparent_resistivity, period, scale='km')[source]: Estimate skin depth from apparent resistivity and period. :param apparent_resistivity: DESCRIPTION. :type apparent_resistivity: TYPE :param period: DESCRIPTION. :type period: TYPE :param scale: DESCRIPTION, defaults to “km”. :type scale: TYPE, optional :return: DESCRIPTION. :rtype: TYPE

from_gocad_sgrid(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

from_modem(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

from_ws3dinv(model_fn)[source]: Read WS3DINV iteration model file. :param model_fn: DESCRIPTION. :type model_fn: TYPE :return: DESCRIPTION. :rtype: TYPE

from_ws3dinv_initial(initial_fn)[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.

Arguments:

**initial_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

get_lower_left_corner(pad_east, pad_north, shift_east=0, shift_north=0)[source]

Get the lower left corner in UTM coordinates for raster. :param shift_north:

Defaults to 0.

Parameters:

shift_east – Defaults to 0.
pad_east (integer) – Number of padding cells to skip from outside in.
pad_north (integer) – Number of padding cells to skip from outside in.

Returns:

Lower left hand corner.

Return type:

mtpy.core.MTLocation

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’

interpolate_to_even_grid(cell_size, pad_north=None, pad_east=None)[source]: Interpolate the model onto an even grid for plotting as a raster or netCDF. :param cell_size: DESCRIPTION. :type cell_size: TYPE :param pad_north: DESCRIPTION, defaults to None. :type pad_north: TYPE, optional :param pad_east: DESCRIPTION, defaults to None. :type pad_east: TYPE, optional :return: DESCRIPTION. :rtype: TYPE

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.

Uses cell_size_east and cell_size_north for dimensions

Adding pad_num of cell_width cells outside of station area

Adding padding cells to given extension and number of padding cells. - extent1 - stretch to a given distance with pad_east or

pad_north number of cells.

extent2 - stretch to a given distance with pad_east or

pad_north number of cells.

stretch stretches from station area using

pad_north and pad_east times pad_stretch_h

Making vertical cells starting with z1_layer increasing logarithmically (base 10) to z_target_depth and num_layers. - default creates a vertical mesh that increases

logarithmically down. See make_z_mesh.

custom input your own vertical mesh.

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.

property save_path: Save path.

to_conductance_raster(cell_size, conductance_dict, pad_north=None, pad_east=None, save_path=None, rotation_angle=0, shift_north=0, shift_east=0, log10=True, verbose=True)[source]

Write out a raster in UTM coordinates for conductance sections.

Expecting a grid that is interoplated onto a regular grid of square cells with size cell_size.

`conductance_dict = {“layer_name”: (min_depth, max_depth)}

if “layer_name” is “surface” then topography is included :param verbose:

Defaults to True.

Parameters:

log10 – Defaults to True.
shift_east – Defaults to 0.
shift_north – Defaults to 0.
cell_size (float) – Square cell size (cell_size x cell_size) in meters.
conductance_dict (TYPE) – DESCRIPTION.
pad_north (integer, optional) – Number of padding cells to skip from outside in, if None defaults to self.pad_north, defaults to None.
pad_east (integer, optional) – Number of padding cells to skip from outside in if None defaults to self.pad_east, defaults to None.
save_path (string or Path, optional) – Path to save files to. If None use self.save_path,, defaults to None.
depth_min (float, optional) – Minimum depth to make raster for in meters,, defaults to None which will use shallowest depth.
depth_max (float, optional) – Maximum depth to make raster for in meters,, defaults to None which will use deepest depth.
rotation_angle (float, optional) – Angle (degrees) to rotate the raster assuming clockwise positive rotation where North = 0, East = 90, defaults to 0.

Raises:

ValueError – If utm_epsg is not input.

Returns:

List of file paths to rasters.

Return type:

TYPE

to_geosoft_xyz(save_fn, 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 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

to_gocad_sgrid(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.

to_modem(model_fn=None, **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:

**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_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'

to_netcdf(fn, pad_east=None, pad_north=None, metadata={})[source]

Create a netCDF file to read into GIS software

works about 50% of the time..

to_raster(cell_size, pad_north=None, pad_east=None, save_path=None, depth_min=None, depth_max=None, rotation_angle=0, shift_north=0, shift_east=0, log10=True, verbose=True)[source]

Write out each depth slice as a raster in UTM coordinates.

Expecting

a grid that is interoplated onto a regular grid of square cells with

size cell_size. :param verbose:

Defaults to True.

param log10:: Defaults to True.
param shift_east:: Defaults to 0.
param shift_north:: Defaults to 0.
param cell_size:: Square cell size (cell_size x cell_size) in meters.
type cell_size:: float
param pad_north:: Number of padding cells to skip from outside in, if None defaults to self.pad_north, defaults to None.
type pad_north:: integer, optional
param pad_east:: Number of padding cells to skip from outside in if None defaults to self.pad_east, defaults to None.
type pad_east:: integer, optional
param save_path:: Path to save files to. If None use self.save_path,, defaults to None.
type save_path:: string or Path, optional
param depth_min:: Minimum depth to make raster for in meters,, defaults to None.
type depth_min:: float, optional
param depth_max:: Maximum depth to make raster for in meters,, defaults to None.
type depth_max:: float, optional
param rotation_angle:: Angle (degrees) to rotate the raster assuming clockwise positive rotation where North = 0, East = 90, defaults to 0.
type rotation_angle:: float, optional
raises ValueError:: If utm_epsg is not input.
return:: List of file paths to rasters.
rtype:: TYPE

to_ubc(basename)[source]: Write a UBC .msh and .mod file. :param basename: :param save_fn: DESCRIPTION. :type save_fn: TYPE :return: DESCRIPTION. :rtype: TYPE

to_vtk(vtk_fn=None, vtk_save_path=None, vtk_fn_basename='ModEM_model_res', **kwargs)[source]

Write a VTK file to plot in 3D rendering programs like Paraview. :param **kwargs: :param vtk_fn: Full path to VKT file to be written, defaults to None. :type vtk_fn: string or Path, optional :param vtk_save_path: Directory to save vtk file to, defaults to None. :type vtk_save_path: string or Path, optional :param vtk_fn_basename: Filename basename of vtk file, note that .vtr

extension is automatically added, defaults to “ModEM_model_res”.

Parameters:

geographic_coordinates (boolean, optional) – [ True | False ] True for geographic coordinates.
units (string, optional) – Units of the spatial grid [ km | m | ft ], defaults to “km”.
coordinate_system – Coordinate system for the station, either the normal MT right-hand coordinate system with z+ down or the sinister z- down [ nez+ | enz- ], defaults to nez+.

Type:

string

Returns:

Full path to VTK file.

Return type:

Path

to_winglink_out(save_fn)[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 :return: DESCRIPTION. :rtype: TYPE

to_ws3dinv_intial(initial_fn, res_list=None)[source]: Write a WS3DINV inital model file.

to_xarray(**kwargs)[source]: Put model in xarray format.

to_xyres(save_path=None, location_type='EN', depth_index='all', outfile_basename='DepthSlice', log_res=False, pad_east=None, pad_north=None)[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.

to_xyzres(savefile=None, location_type='EN', log_res=False, pad_east=None, pad_north=None)[source]: Save a model file as a space delimited x y z res file.

mtpy.modeling.winglink module

Created on Mon Aug 22 15:19:30 2011

deal with output files from winglink.

@author: jp

class mtpy.modeling.winglink.PlotMisfitPseudoSection(data_fn, resp_fn, **kwargs)[source]

Bases: object

plot a pseudo section misfit of the data and response if given

Note

the output file from winglink does not contain errors, so to get a normalized error, you need to input the error for each component as a percent for resistivity and a value for phase and tipper. If you used the data errors, unfortunately, you have to input those as arrays.

wl_data_fnstring
full path to output data file from winglink

key words

description

axmpte

matplotlib.axes instance for TE model phase

axmptm

matplotlib.axes instance for TM model phase

axmrte

matplotlib.axes instance for TE model app. res

axmrtm

matplotlib.axes instance for TM model app. res

axpte

matplotlib.axes instance for TE data phase

axptm

matplotlib.axes instance for TM data phase

axrte

matplotlib.axes instance for TE data app. res.

axrtm

matplotlib.axes instance for TM data app. res.

cb_pad

padding between colorbar and axes

cb_shrink

percentage to shrink the colorbar to

fig

matplotlib.figure instance

fig_dpi

resolution of figure in dots per inch

fig_num

number of figure instance

fig_size

size of figure in inches (width, height)

font_size

size of font in points

label_list

list to label plots

ml

factor to label stations if 2 every other station is labeled on the x-axis

period

np.array of periods to plot

phase_cmap

color map name of phase

phase_limits_te

limits for te phase in degrees (min, max)

phase_limits_tm

limits for tm phase in degrees (min, max)

plot_resp

[ ‘y’ | ‘n’ ] to plot response

plot_yn

[ ‘y’ | ‘n’ ] ‘y’ to plot on instantiation

res_cmap

color map name for resistivity

res_limits_te

limits for te resistivity in log scale (min, max)

res_limits_tm

limits for tm resistivity in log scale (min, max)

rp_list

list of dictionaries as made from read2Dresp

station_id

index to get station name (min, max)

station_list

station list got from rp_list

subplot_bottom

subplot spacing from bottom (relative coordinates)

subplot_hspace

vertical spacing between subplots

subplot_left

subplot spacing from left

subplot_right

subplot spacing from right

subplot_top

subplot spacing from top

subplot_wspace

horizontal spacing between subplots

Methods

Description

plot

plots a pseudo-section of apparent resistiviy and phase of data and model if given. called on instantiation if plot_yn is ‘y’.

redraw_plot

call redraw_plot to redraw the figures, if one of the attributes has been changed

save_figure

saves the matplotlib.figure instance to desired location and format

Example:

>>> import mtpy.modeling.occam2d as occam2d
>>> ocd = occam2d.Occam2DData()
>>> rfile = r"/home/Occam2D/Line1/Inv1/Test_15.resp"
>>> ocd.data_fn = r"/home/Occam2D/Line1/Inv1/DataRW.dat"
>>> ps1 = ocd.plot2PseudoSection(resp_fn=rfile)

get_misfit()[source]

compute misfit of MT response found from the model and the data.

Need to normalize correctly

plot()[source]: plot pseudo section of data and response if given

redraw_plot()[source]

redraw plot if parameters were changed

use this function if you updated some attributes and want to re-plot.

Example:

>>> # change the color and marker of the xy components
>>> import mtpy.modeling.occam2d as occam2d
>>> ocd = occam2d.Occam2DData(r"/home/occam2d/Data.dat")
>>> p1 = ocd.plotPseudoSection()
>>> #change color of te markers to a gray-blue
>>> p1.res_cmap = 'seismic_r'
>>> p1.redraw_plot()

save_figure(save_fn, file_format='pdf', orientation='portrait', fig_dpi=None, close_plot='y')[source]

save_plot will save the figure to save_fn.

Arguments:

save_fnstring
full path to save figure to, can be input as * directory path -> the directory path to save to

in which the file will be saved as save_fn/station_name_PhaseTensor.file_format

full path -> file will be save to the given path. If you use this option then the format will be assumed to be provided by the path

file_format[ pdf | eps | jpg | png | svg ]
file type of saved figure pdf,svg,eps…

orientation[ landscape | portrait ]
orientation in which the file will be saved default is portrait

fig_dpiint
The resolution in dots-per-inch the file will be saved. If None then the dpi will be that at which the figure was made. I don’t think that it can be larger than dpi of the figure.

close_plot[ y | n ]

‘y’ will close the plot after saving.

‘n’ will leave plot open

Example:

>>> # to save plot as jpg
>>> import mtpy.modeling.occam2d as occam2d
>>> dfn = r"/home/occam2d/Inv1/data.dat"
>>> ocd = occam2d.Occam2DData(dfn)
>>> ps1 = ocd.plotPseudoSection()
>>> ps1.save_plot(r'/home/MT/figures', file_format='jpg')

update_plot()[source]

update any parameters that where changed using the built-in draw from canvas.

Use this if you change an of the .fig or axes properties

Example:

>>> # to change the grid lines to only be on the major ticks
>>> import mtpy.modeling.occam2d as occam2d
>>> dfn = r"/home/occam2d/Inv1/data.dat"
>>> ocd = occam2d.Occam2DData(dfn)
>>> ps1 = ocd.plotPseudoSection()
>>> [ax.grid(True, which='major') for ax in [ps1.axrte,ps1.axtep]]
>>> ps1.update_plot()

class mtpy.modeling.winglink.PlotPseudoSection(wl_data_fn=None, **kwargs)[source]

Bases: object

plot a pseudo section of the data and response if given

wl_data_fnstring
full path to winglink output data file.

key words

description

axmpte

matplotlib.axes instance for TE model phase

axmptm

matplotlib.axes instance for TM model phase

axmrte

matplotlib.axes instance for TE model app. res

axmrtm

matplotlib.axes instance for TM model app. res

axpte

matplotlib.axes instance for TE data phase

axptm

matplotlib.axes instance for TM data phase

axrte

matplotlib.axes instance for TE data app. res.

axrtm

matplotlib.axes instance for TM data app. res.

cb_pad

padding between colorbar and axes

cb_shrink

percentage to shrink the colorbar to

fig

matplotlib.figure instance

fig_dpi

resolution of figure in dots per inch

fig_num

number of figure instance

fig_size

size of figure in inches (width, height)

font_size

size of font in points

label_list

list to label plots

ml

factor to label stations if 2 every other station is labeled on the x-axis

period

np.array of periods to plot

phase_cmap

color map name of phase

phase_limits_te

limits for te phase in degrees (min, max)

phase_limits_tm

limits for tm phase in degrees (min, max)

plot_resp

[ ‘y’ | ‘n’ ] to plot response

plot_tipper

[ ‘y’ | ‘n’ ] to plot tipper

plot_yn

[ ‘y’ | ‘n’ ] ‘y’ to plot on instantiation

res_cmap

color map name for resistivity

res_limits_te

limits for te resistivity in log scale (min, max)

res_limits_tm

limits for tm resistivity in log scale (min, max)

rp_list

list of dictionaries as made from read2Dresp

station_id

index to get station name (min, max)

station_list

station list got from rp_list

subplot_bottom

subplot spacing from bottom (relative coordinates)

subplot_hspace

vertical spacing between subplots

subplot_left

subplot spacing from left

subplot_right

subplot spacing from right

subplot_top

subplot spacing from top

subplot_wspace

horizontal spacing between subplots

Methods

Description

plot

plots a pseudo-section of apparent resistiviy and phase of data and model if given. called on instantiation if plot_yn is ‘y’.

redraw_plot

call redraw_plot to redraw the figures, if one of the attributes has been changed

save_figure

saves the matplotlib.figure instance to desired location and format

Example:

>>> import mtpy.modeling.winglink as winglink
>>> d_fn = r"/home/winglink/Line1/Inv1/DataRW.txt"
>>> ps_plot = winglink.PlotPseudoSection(d_fn)

plot()[source]: plot pseudo section of data and response if given

redraw_plot()[source]

redraw plot if parameters were changed

use this function if you updated some attributes and want to re-plot.

Example:

>>> # plot tipper and change station id
>>> import mtpy.modeling.winglink as winglink
>>> ps_plot = winglink.PlotPseudosection(wl_fn)
>>> ps_plot.plot_tipper = 'y'
>>> ps_plot.station_id = [2, 5]
>>> #label only every 3rd station
>>> ps_plot.ml = 3
>>> ps_plot.redraw_plot()

save_figure(save_fn, file_format='pdf', orientation='portrait', fig_dpi=None, close_plot='y')[source]

save_plot will save the figure to save_fn.

Arguments:

save_fnstring
full path to save figure to, can be input as * directory path -> the directory path to save to

in which the file will be saved as save_fn/station_name_PhaseTensor.file_format

full path -> file will be save to the given path. If you use this option then the format will be assumed to be provided by the path

file_format[ pdf | eps | jpg | png | svg ]
file type of saved figure pdf,svg,eps…

orientation[ landscape | portrait ]
orientation in which the file will be saved default is portrait

fig_dpiint
The resolution in dots-per-inch the file will be saved. If None then the dpi will be that at which the figure was made. I don’t think that it can be larger than dpi of the figure.

close_plot[ y | n ]

‘y’ will close the plot after saving.

‘n’ will leave plot open

Example:

>>> # to save plot as jpg
>>> ps_plot.save_plot(r'/home/MT/figures', file_format='jpg')

update_plot()[source]

update any parameters that where changed using the built-in draw from canvas.

Use this if you change an of the .fig or axes properties

Example:

>>> # to change the grid lines to only be on the major ticks
>>> [ax.grid(True, which='major') for ax in [ps_plot.axrte]]
>>> ps_plot.update_plot()

class mtpy.modeling.winglink.PlotResponse(wl_data_fn=None, resp_fn=None, **kwargs)[source]

Bases: object

Helper class to deal with plotting the MT response and occam2d model.

Arguments:

data_fnstring
full path to data file

resp_fnstring or list
full path(s) to response file(s)

Attributes/key words	description
ax_list	list of matplotlib.axes instances for use with OccamPointPicker
color_mode	[ ‘color’ \| ‘bw’ ] plot figures in color or black and white (‘bw’)
cted	color of Data TE marker and line
ctem	color of Model TE marker and line
ctewl	color of Winglink Model TE marker and line
ctmd	color of Data TM marker and line
ctmm	color of Model TM marker and line
ctmwl	color of Winglink Model TM marker and line
e_capsize	size of error bar caps in points
e_capthick	line thickness of error bar caps in points
err_list	list of line properties of error bars for use with OccamPointPicker
fig_dpi	figure resolution in dots-per-inch
fig_list	list of dictionaries with key words station –> station name fig –> matplotlib.figure instance axrte –> matplotlib.axes instance for TE app.res axrtm –> matplotlib.axes instance for TM app.res axpte –> matplotlib.axes instance for TE phase axptm –> matplotlib.axes instance for TM phase
fig_num	starting number of figure
fig_size	size of figure in inches (width, height)
font_size	size of axes ticklabel font in points
line_list	list of matplotlib.Line instances for use with OccamPointPicker
lw	line width of lines in points
ms	marker size in points
mted	marker for Data TE mode
mtem	marker for Model TE mode
mtewl	marker for Winglink Model TE
mtmd	marker for Data TM mode
mtmm	marker for Model TM mode
mtmwl	marker for Winglink TM mode
period	np.ndarray of periods to plot
phase_limits	limits on phase plots in degrees (min, max)
plot_model_error	[ ‘y’ \| ‘n’ ] default is ‘y’ to plot model errors
plot_num	[ 1 \| 2 ] 1 to plot both modes in a single plot 2 to plot modes in separate plots (default)
plot_tipper	[ ‘y’ \| ‘n’ ] plot tipper data if desired
plot_type	[ ‘1’ \| station_list] ‘1’ –> to plot all stations in different figures station_list –> to plot a few stations, give names of stations ex. [‘mt01’, ‘mt07’]
plot_yn	[ ‘y’ \| ‘n’] ‘y’ –> to plot on instantiation ‘n’ –> to not plot on instantiation
res_limits	limits on resistivity plot in log scale (min, max)
rp_list	list of dictionaries from read2Ddata
station_list	station_list list of stations in rp_list
subplot_bottom	subplot spacing from bottom (relative coordinates)
subplot_hspace	vertical spacing between subplots
subplot_left	subplot spacing from left
subplot_right	subplot spacing from right
subplot_top	subplot spacing from top
subplot_wspace	horizontal spacing between subplots
wl_fn	Winglink file name (full path)

Methods	Description
plot	plots the apparent resistiviy and phase of data and model if given. called on instantiation if plot_yn is ‘y’.
redraw_plot	call redraw_plot to redraw the figures, if one of the attributes has been changed
save_figures	save all the matplotlib.figure instances in fig_list

Example:

:: >>> data_fn = r”/home/occam/line1/inv1/OccamDataFile.dat” >>> resp_list = [r”/home/occam/line1/inv1/test_{0:02}”.format(ii)

for ii in range(2, 8, 2)]

>>> pr_obj = occam2d.PlotResponse(data_fn, resp_list, plot_tipper='y')

plot()[source]: plot the data and model response, if given, in individual plots.

redraw_plot()[source]

redraw plot if parameters were changed

use this function if you updated some attributes and want to re-plot.

Example:

>>> # change the color and marker of the xy components
>>> import mtpy.modeling.occam2d as occam2d
>>> ocd = occam2d.Occam2DData(r"/home/occam2d/Data.dat")
>>> p1 = ocd.plot2DResponses()
>>> #change color of te markers to a gray-blue
>>> p1.cted = (.5, .5, .7)
>>> p1.redraw_plot()

save_figures(save_path, fig_fmt='pdf', fig_dpi=None, close_fig='y')[source]

save all the figure that are in self.fig_list

Example:

>>> # change the color and marker of the xy components
>>> import mtpy.modeling.occam2d as occam2d
>>> ocd = occam2d.Occam2DData(r"/home/occam2d/Data.dat")
>>> p1 = ocd.plot2DResponses()
>>> p1.save_figures(r"/home/occam2d/Figures", fig_fmt='jpg')

exception mtpy.modeling.winglink.WLInputError[source]: Bases: Exception

mtpy.modeling.winglink.read_model_file(model_fn)[source]

readModelFile reads in the XYZ txt file output by Winglink.

Inputs:: modelfile = fullpath and filename to modelfile profiledirection = ‘ew’ for east-west predominantly, ‘ns’ for

predominantly north-south. This gives column to fix

mtpy.modeling.winglink.read_output_file(output_fn)[source]

Reads in an output file from winglink and returns the data in the form of a dictionary of structured arrays.

Arguments:

output_fnstring
the full path to winglink outputfile

Returns:

wl_datadictionary with keys of station names
each station contains a structured array with keys * ‘station’ –> station name * ‘period’ –> periods to plot * ‘te_res’ –> TE resistivity in linear scale * ‘tm_res’ –> TM resistivity in linear scale * ‘te_phase’ –> TE phase in deg * ‘tm_phase’ –> TM phase in deg * ‘re_tip’ –> real tipper amplitude. * ‘im_tip’ –> imaginary tipper amplitude * ‘rms’ –> RMS for the station * ‘index’ –> order from left to right of station number

Note

each data is an np.ndarray(2, num_periods) where the first index is the data and the second index is the model response

mtpy.modeling.winglinktools module

Created on Mon Aug 22 15:19:30 2011

@author: a1185872

mtpy.modeling.winglinktools.plotResponses(outputfile, maxcol=8, plottype='all', **kwargs)[source]

PlotResponse will plot the responses modeled from winglink against the observed data.

Inputs:: outputfile = full path and filename to output file maxcol = maximum number of columns for the plot plottype = ‘all’ to plot all on the same plot

‘1’ to plot each respones in a different figure station to plot a single station or enter as a list of stations to plot a few stations [station1,station2]. Does not have to be verbatim but should have similar unique characters input pb01 for pb01cs in outputfile
Outputs:: None

mtpy.modeling.winglinktools.readModelFile(modelfile, profiledirection='ew')[source]

ReadModelFile reads in the XYZ txt file output by Winglink.

Inputs:: modelfile = fullpath and filename to modelfile profiledirection = ‘ew’ for east-west predominantly, ‘ns’ for

predominantly north-south. This gives column to fix.

mtpy.modeling.winglinktools.readOutputFile(outputfile)[source]

ReadOutputFile will read an output file from winglink and output data in the form of a dictionary.

Input:

outputfile = the full path and filename of outputfile

Output:

idict = dictionary with keys of station name: each idict[station name] is a dictionary with keys corresponding to modeled and observed responses:

‘obsresxy’,’obsphasexy’,’modresxy’,’modphasexy’,’obsresyx’, ‘obsphaseyx’,’modresyx’,’modphaseyx’,’obshzres’, ‘obshzphase’,’modhzres’,’modhzphase’,’period’
rplst = list of dictionaries for each station with keywords:: ‘station’ = station name ‘offset’ = relative offset, ‘resxy’ = TE resistivity and error as row 0 and 1 respectively, ‘resyx’= TM resistivity and error as row 0 and 1 respectively, ‘phasexy’= TE phase and error as row 0 and 1 respectively, ‘phaseyx’= Tm phase and error as row 0 and 1 respectively, ‘realtip’= Real Tipper and error as row 0 and 1 respectively, ‘imagtip’= Imaginary Tipper and error as row 0 and 1 respectively

plst = periodlst as the median of all stations. stationlst = list of stations in order from profile title = list of parameters for plotting as [title,profile,inversiontype]

Module contents

class mtpy.modeling.StructuredGrid3D(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

convert_model_to_int(res_list=None)[source]: Convert resistivity values to integers according to resistivity list. :param res_list: Resistivity values in Ohm-m, defaults to None. :type res_list: list of floats, optional :return: Array of integers corresponding to the res_list. :rtype: np.ndarray(dtype=int)

estimate_skin_depth(apparent_resistivity, period, scale='km')[source]: Estimate skin depth from apparent resistivity and period. :param apparent_resistivity: DESCRIPTION. :type apparent_resistivity: TYPE :param period: DESCRIPTION. :type period: TYPE :param scale: DESCRIPTION, defaults to “km”. :type scale: TYPE, optional :return: DESCRIPTION. :rtype: TYPE

from_gocad_sgrid(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

from_modem(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

from_ws3dinv(model_fn)[source]: Read WS3DINV iteration model file. :param model_fn: DESCRIPTION. :type model_fn: TYPE :return: DESCRIPTION. :rtype: TYPE

from_ws3dinv_initial(initial_fn)[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.

Arguments:

**initial_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

get_lower_left_corner(pad_east, pad_north, shift_east=0, shift_north=0)[source]

Get the lower left corner in UTM coordinates for raster. :param shift_north:

Defaults to 0.

Parameters:

shift_east – Defaults to 0.
pad_east (integer) – Number of padding cells to skip from outside in.
pad_north (integer) – Number of padding cells to skip from outside in.