Package Modeling¶
Module ModEM¶
- class mtpy.modeling.modem.ControlFwd(**kwargs)[source]¶
read and write control file for
This file controls how the inversion starts and how it is run
Methods
read_control_file
([control_fn])read in a control file
write_control_file
([control_fn, save_path, ...])write control file
- class mtpy.modeling.modem.ControlInv(**kwargs)[source]¶
read and write control file for how the inversion starts and how it is run
Methods
read_control_file
([control_fn])read in a control file
write_control_file
([control_fn, save_path, ...])write control file
- class mtpy.modeling.modem.Covariance(grid_dimensions=None, **kwargs)[source]¶
read and write covariance files
Methods
read_cov_file
(cov_fn)read a covariance file
write_cov_vtk_file
(cov_vtk_fn[, model_fn, ...])write a vtk file of the covariance to match things up
write_covariance_file
([cov_fn, save_path, ...])write a covariance file
get_parameters
- class mtpy.modeling.modem.Data(edi_list=None, **kwargs)[source]¶
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
- Attributes
rotation_angle
Rotate data assuming N=0, E=90
station_locations
location of stations
Methods
center_stations
(model_fn[, data_fn])Center station locations to the middle of cells, might be useful for topography.
change_data_elevation
(model_obj[, data_fn, ...])At each station in the data file rewrite the elevation, so the station is on the surface, not floating in air.
compute the error from the given parameters
compute_phase_tensor
(datfile, outdir)Compute the phase tensors from a ModEM dat file :param datfile: path2/file.dat :return: path2csv created by this method
convert_modem_to_ws
([data_fn, ws_data_fn, ...])convert a ModEM data file to WS format.
convert_ws3dinv_data_file
(ws_data_fn[, ...])convert a ws3dinv data file into ModEM format
fill_data_array
([new_edi_dir, ...])fill the data array from mt_dict
filter_periods
(mt_obj, per_array)Select the periods of the mt_obj that are in per_array.
get_header_string
(error_type, error_value, ...)reset the header sring for file
get mt_dict from edi file list
get important parameters for documentation
make a period list to invert for
get station locations from edi files
project_stations_on_topography
(model_object)Re-write the data file to change the elevation column.
read_data_file
([data_fn, center_utm])Read ModEM data file
write_data_file
([save_path, fn_basename, ...])write data file for ModEM will save file as save_path/fn_basename
write_vtk_station_file
([vtk_save_path, ...])write a vtk file for station locations.
- center_stations(model_fn, data_fn=None)[source]¶
Center station locations to the middle of cells, might be useful for topography.
- Returns
- **new_data_fn**string
full path to new data file
- change_data_elevation(model_obj, data_fn=None, res_air=1000000000000.0)[source]¶
At each station in the data file rewrite the elevation, so the station is on the surface, not floating in air.
- compute_phase_tensor(datfile, outdir)[source]¶
Compute the phase tensors from a ModEM dat file :param datfile: path2/file.dat :return: path2csv created by this method
- convert_modem_to_ws(data_fn=None, ws_data_fn=None, error_map=[1, 1, 1, 1])[source]¶
convert a ModEM data file to WS format.
- convert_ws3dinv_data_file(ws_data_fn, station_fn=None, save_path=None, fn_basename=None)[source]¶
convert a ws3dinv data file into ModEM format
- fill_data_array(new_edi_dir=None, use_original_freq=False, longitude_format='LON')[source]¶
fill the data array from mt_dict
- static filter_periods(mt_obj, per_array)[source]¶
Select the periods of the mt_obj that are in per_array. used to do original freq inversion.
- Parameters
mt_obj –
per_array –
- Returns
array of selected periods (subset) of the mt_obj
- static get_header_string(error_type, error_value, rotation_angle)[source]¶
reset the header sring for file
- project_stations_on_topography(model_object, air_resistivity=1000000000000.0)[source]¶
Re-write the data file to change the elevation column. And update covariance mask according topo elevation model. :param model_object: :param air_resistivity: :return:
- read_data_file(data_fn=None, center_utm=None)[source]¶
Read ModEM data file
- inputs:
data_fn = full path to data file name center_utm = option to provide real world coordinates of the center of
the grid for putting the data and model back into utm/grid coordinates, format [east_0, north_0, z_0]
- Fills attributes:
data_array
period_list
mt_dict
- property rotation_angle¶
Rotate data assuming N=0, E=90
- property station_locations¶
location of stations
- class mtpy.modeling.modem.ModEMConfig(**kwargs)[source]¶
read and write configuration files for how each inversion is run
Methods
add_dict
([fn, obj])add dictionary based on file name or object
write_config_file
([save_dir, config_fn_basename])write a config file based on provided information
- class mtpy.modeling.modem.Model(stations_object=None, data_object=None, **kwargs)[source]¶
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.
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
- Attributes
- nodes_east
- nodes_north
- nodes_z
- plot_east
- plot_north
- plot_z
Methods
add_layers_to_mesh
([n_add_layers, ...])Function to add constant thickness layers to the top or bottom of mesh.
add_topography_from_data
(data_object[, ...])Wrapper around add_topography_to_model2 that allows creating a surface model from EDI data.
add_topography_to_model2
([topographyfile, ...])if air_layers is non-zero, will add topo: read in topograph file, make a surface model.
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)
get important model parameters to write to a file for documentation later.
interpolate_elevation2
([surfacefile, ...])project a surface to the model grid and add resulting elevation data to a dictionary called surface_dict.
create finite element mesh according to user-input parameters.
make_z_mesh_new
([n_layers])new version of make_z_mesh.
plot_mesh
([east_limits, north_limits, z_limits])Plot the mesh to show model grid
# add mesh grid lines in xy plan north-east map :return:
display the mesh in North-Depth aspect :return:
create a quick pcolor plot of the resistivity at sea level with stations, to check if we have stations in the sea
display topography elevation data together with station locations on a cell-index N-E map :return:
read_gocad_sgrid_file
(sgrid_header_file[, ...])read a gocad sgrid file and put this info into a ModEM file.
read_model_file
([model_fn])read an initial file and return the pertinent information including grid positions in coordinates relative to the center point (0,0) and starting model.
read_ws_model_file
(ws_model_fn)reads in a WS3INV3D model file
write_gocad_sgrid_file
([fn, origin, clip, ...])write a model to gocad sgrid
write_model_file
(**kwargs)will write an initial file for ModEM.
write_vtk_file
([vtk_save_path, vtk_fn_basename])write a vtk file to view in Paraview or other
write_xyres
([savepath, location_type, ...])write files containing depth slice data (x, y, res for each depth)
write_xyzres
([savefile, location_type, ...])save a model file as a space delimited x y z res file
print_mesh_params
print_model_file_summary
- 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_model2
- Parameters
n_add_layers – integer, number of layers to add
layer_thickness – real value or list/array. Thickness of layers, defaults to z1 layer. Can provide a single value or a list/array containing multiple layer thicknesses.
where – where to add, top or bottom
- add_topography_from_data(data_object, interp_method='nearest', air_resistivity=1000000000000.0, topography_buffer=None, airlayer_type='log_up')[source]¶
Wrapper around add_topography_to_model2 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_model2 so a surface model can be interpolated from it.
The surface tuple is of format (lon, lat, elev) containing station locations.
- Args:
- data_object (mtpy.modeling.ModEM.data.Data): A ModEm data
object that has been filled with data from EDI files.
- interp_method (str, optional): Same as
add_topography_to_model2.
- air_resistivity (float, optional): Same as
add_topography_to_model2.
- topography_buffer (float): Same as
add_topography_to_model2.
- airlayer_type (str, optional): Same as
add_topography_to_model2.
- add_topography_to_model2(topographyfile=None, surface=None, topographyarray=None, interp_method='nearest', air_resistivity=1000000000000.0, topography_buffer=None, airlayer_type='log_up', max_elev=None)[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.
- Parameters
topographyfile – file containing topography (arcgis ascii grid)
topographyarray – alternative to topographyfile - array of elevation values on model grid
interp_method – interpolation method for topography, ‘nearest’, ‘linear’, or ‘cubic’
air_resistivity – resistivity value to assign to air
topography_buffer – buffer around stations to calculate minimum and maximum topography value to use for meshing
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
- assign_resistivity_from_surfacedata(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 surfacename = 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
- get_parameters()[source]¶
get important model parameters to write to a file for documentation later.
- interpolate_elevation2(surfacefile=None, surface=None, get_surfacename=False, method='nearest', fast=True)[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 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: 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()[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
- plot_mesh(east_limits=None, north_limits=None, z_limits=None, **kwargs)[source]¶
Plot the mesh to show model grid
- plot_sealevel_resistivity()[source]¶
create a quick pcolor plot of the resistivity at sea level with stations, to check if we have stations in the sea
- plot_topography()[source]¶
display topography elevation data together with station locations on a cell-index N-E map :return:
- 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
- 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.
- write_vtk_file(vtk_save_path=None, vtk_fn_basename='ModEM_model_res')[source]¶
write a vtk file to view in Paraview or other
- write_xyres(savepath=None, location_type='EN', origin=[0, 0], model_epsg=None, depth_index='all', outfile_basename='DepthSlice', log_res=False, model_utm_zone=None, 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)
savepath = 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
- class mtpy.modeling.modem.ModelManipulator(model_fn=None, data_fn=None, **kwargs)[source]¶
will plot a model from wsinv3d or init file so the user can manipulate the resistivity values relatively easily. At the moment only plotted in map view.
- Example
:: >>> import mtpy.modeling.ws3dinv as ws >>> initial_fn = r”/home/MT/ws3dinv/Inv1/WSInitialFile” >>> mm = ws.WSModelManipulator(initial_fn=initial_fn)
Buttons
Description
‘=’
increase depth to next vertical node (deeper)
‘-’
decrease depth to next vertical node (shallower)
‘q’
quit the plot, rewrites initial file when pressed
‘a’
copies the above horizontal layer to the present layer
‘b’
copies the below horizonal layer to present layer
‘u’
undo previous change
Attributes
Description
ax1
matplotlib.axes instance for mesh plot of the model
ax2
matplotlib.axes instance of colorbar
cb
matplotlib.colorbar instance for colorbar
cid_depth
matplotlib.canvas.connect for depth
cmap
matplotlib.colormap instance
cmax
maximum value of resistivity for colorbar. (linear)
cmin
minimum value of resistivity for colorbar (linear)
data_fn
full path fo data file
depth_index
integer value of depth slice for plotting
dpi
resolution of figure in dots-per-inch
dscale
depth scaling, computed internally
east_line_xlist
list of east mesh lines for faster plotting
east_line_ylist
list of east mesh lines for faster plotting
fdict
dictionary of font properties
fig
matplotlib.figure instance
fig_num
number of figure instance
fig_size
size of figure in inches
font_size
size of font in points
grid_east
location of east nodes in relative coordinates
grid_north
location of north nodes in relative coordinates
grid_z
location of vertical nodes in relative coordinates
initial_fn
full path to initial file
m_height
mean height of horizontal cells
m_width
mean width of horizontal cells
map_scale
[ ‘m’ | ‘km’ ] scale of map
mesh_east
np.meshgrid of east, north
mesh_north
np.meshgrid of east, north
mesh_plot
matplotlib.axes.pcolormesh instance
model_fn
full path to model file
new_initial_fn
full path to new initial file
nodes_east
spacing between east nodes
nodes_north
spacing between north nodes
nodes_z
spacing between vertical nodes
north_line_xlist
list of coordinates of north nodes for faster plotting
north_line_ylist
list of coordinates of north nodes for faster plotting
plot_yn
[ ‘y’ | ‘n’ ] plot on instantiation
radio_res
matplotlib.widget.radio instance for change resistivity
rect_selector
matplotlib.widget.rect_selector
res
np.ndarray(nx, ny, nz) for model in linear resistivity
res_copy
copy of res for undo
res_dict
dictionary of segmented resistivity values
res_list
list of resistivity values for model linear scale
res_model
np.ndarray(nx, ny, nz) of resistivity values from res_list (linear scale)
res_model_int
np.ndarray(nx, ny, nz) of integer values corresponding to res_list for initial model
res_value
current resistivty value of radio_res
save_path
path to save initial file to
station_east
station locations in east direction
station_north
station locations in north direction
xlimits
limits of plot in e-w direction
ylimits
limits of plot in n-s direction
- Attributes
- nodes_east
- nodes_north
- nodes_z
- plot_east
- plot_north
- plot_z
Methods
add_layers_to_mesh
([n_add_layers, ...])Function to add constant thickness layers to the top or bottom of mesh.
add_topography_from_data
(data_object[, ...])Wrapper around add_topography_to_model2 that allows creating a surface model from EDI data.
add_topography_to_model2
([topographyfile, ...])if air_layers is non-zero, will add topo: read in topograph file, make a surface model.
assign_resistivity_from_surfacedata
(...)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)
change_model_res
(xchange, ychange)change resistivity values of resistivity model
reads in initial file or model file and set attributes:
get_parameters
()get important model parameters to write to a file for documentation later.
interpolate_elevation2
([surfacefile, ...])project a surface to the model grid and add resulting elevation data to a dictionary called surface_dict.
make_mesh
()create finite element mesh according to user-input parameters.
make_z_mesh_new
([n_layers])new version of make_z_mesh.
plot
()plots the model with:
plot_mesh
([east_limits, north_limits, z_limits])Plot the mesh to show model grid
plot_mesh_xy
()# add mesh grid lines in xy plan north-east map :return:
plot_mesh_xz
()display the mesh in North-Depth aspect :return:
plot_sealevel_resistivity
()create a quick pcolor plot of the resistivity at sea level with stations, to check if we have stations in the sea
plot_topography
()display topography elevation data together with station locations on a cell-index N-E map :return:
read_gocad_sgrid_file
(sgrid_header_file[, ...])read a gocad sgrid file and put this info into a ModEM file.
read_model_file
([model_fn])read an initial file and return the pertinent information including grid positions in coordinates relative to the center point (0,0) and starting model.
read_ws_model_file
(ws_model_fn)reads in a WS3INV3D model file
rect_onselect
(eclick, erelease)on selecting a rectangle change the colors to the resistivity values
redraws the plot
rewrite_model_file
([model_fn, save_path, ...])write an initial file for wsinv3d from the model created.
set_res_list
(res_list)on setting res_list also set the res_dict to correspond
set_res_value
(val)write_gocad_sgrid_file
([fn, origin, clip, ...])write a model to gocad sgrid
write_model_file
(**kwargs)will write an initial file for ModEM.
write_vtk_file
([vtk_save_path, vtk_fn_basename])write a vtk file to view in Paraview or other
write_xyres
([savepath, location_type, ...])write files containing depth slice data (x, y, res for each depth)
write_xyzres
([savefile, location_type, ...])save a model file as a space delimited x y z res file
print_mesh_params
print_model_file_summary
- get_model()[source]¶
- reads in initial file or model file and set attributes:
-resmodel -northrid -eastrid -zgrid -res_list if initial file
- plot()[source]¶
- plots the model with:
-a radio dial for depth slice -radio dial for resistivity value
- rect_onselect(eclick, erelease)[source]¶
on selecting a rectangle change the colors to the resistivity values
- class mtpy.modeling.modem.PlotResponse(data_fn=None, resp_fn=None, **kwargs)[source]¶
plot data and response
Plots the real and imaginary impedance and induction vector if present.
- Example
>>> import mtpy.modeling.modem as modem >>> dfn = r"/home/MT/ModEM/Inv1/DataFile.dat" >>> rfn = r"/home/MT/ModEM/Inv1/Test_resp_000.dat" >>> mrp = modem.PlotResponse(data_fn=dfn, resp_fn=rfn) >>> # plot only the TE and TM modes >>> mrp.plot_component = 2 >>> mrp.redraw_plot()
Attributes
Description
color_mode
[ ‘color’ | ‘bw’ ] color or black and white plots
cted
color for data Z_XX and Z_XY mode
ctem
color for model Z_XX and Z_XY mode
ctmd
color for data Z_YX and Z_YY mode
ctmm
color for model Z_YX and Z_YY mode
data_fn
full path to data file
data_object
WSResponse instance
e_capsize
cap size of error bars in points (default is .5)
e_capthick
cap thickness of error bars in points (default is 1)
fig_dpi
resolution of figure in dots-per-inch (300)
fig_list
list of matplotlib.figure instances for plots
fig_size
size of figure in inches (default is [6, 6])
font_size
size of font for tick labels, axes labels are font_size+2 (default is 7)
legend_border_axes_pad
padding between legend box and axes
legend_border_pad
padding between border of legend and symbols
legend_handle_text_pad
padding between text labels and symbols of legend
legend_label_spacing
padding between labels
legend_loc
location of legend
legend_marker_scale
scale of symbols in legend
lw
line width data curves (default is .5)
ms
size of markers (default is 1.5)
lw_r
line width response curves (default is .5)
ms_r
size of markers response curves (default is 1.5)
mted
marker for data Z_XX and Z_XY mode
mtem
marker for model Z_XX and Z_XY mode
mtmd
marker for data Z_YX and Z_YY mode
mtmm
marker for model Z_YX and Z_YY mode
phase_limits
limits of phase
plot_component
[ 2 | 4 ] 2 for TE and TM or 4 for all components
plot_style
[ 1 | 2 ] 1 to plot each mode in a seperate subplot and 2 to plot xx, xy and yx, yy in same plots
plot_type
[ ‘1’ | list of station name ] ‘1’ to plot all stations in data file or input a list of station names to plot if station_fn is input, otherwise input a list of integers associated with the index with in the data file, ie 2 for 2nd station
plot_z
[ True | False ] default is True to plot impedance, False for plotting resistivity and phase
plot_yn
[ ‘n’ | ‘y’ ] to plot on instantiation
res_limits
limits of resistivity in linear scale
resp_fn
full path to response file
resp_object
WSResponse object for resp_fn, or list of WSResponse objects if resp_fn is a list of response files
station_fn
full path to station file written by WSStation
subplot_bottom
space between axes and bottom of figure
subplot_hspace
space between subplots in vertical direction
subplot_left
space between axes and left of figure
subplot_right
space between axes and right of figure
subplot_top
space between axes and top of figure
subplot_wspace
space between subplots in horizontal direction
Methods
redraw plot if parameters were changed
save_figure
(save_fn[, file_format, ...])save_plot will save the figure to save_fn.
plot
- 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.plotAllResponses() >>> #change line width >>> p1.lw = 2 >>> p1.redraw_plot()
- class mtpy.modeling.modem.PlotSlices(model_fn, data_fn=None, **kwargs)[source]¶
Plot all cartesian axis-aligned slices and be able to scroll through the model
Extract arbitrary profiles (e.g. along a seismic line) from a model
- Example
>>> import mtpy.modeling.modem as modem >>> mfn = r"/home/modem/Inv1/Modular_NLCG_100.rho" >>> dfn = r"/home/modem/Inv1/ModEM_data.dat" >>> pds = ws.PlotSlices(model_fn=mfn, data_fn=dfn)
Buttons
Description
‘e’
moves n-s slice east by one model block
‘w’
moves n-s slice west by one model block
‘n’
moves e-w slice north by one model block
‘m’
moves e-w slice south by one model block
‘d’
moves depth slice down by one model block
‘u’
moves depth slice up by one model block
Attributes
Description
ax_en
matplotlib.axes instance for depth slice map view
ax_ez
matplotlib.axes instance for e-w slice
ax_map
matplotlib.axes instance for location map
ax_nz
matplotlib.axes instance for n-s slice
climits
(min , max) color limits on resistivity in log scale. default is (0, 4)
cmap
name of color map for resisitiviy. default is ‘jet_r’
data_fn
full path to data file name
draw_colorbar
show colorbar on exported plot; default True
dscale
scaling parameter depending on map_scale
east_line_xlist
list of line nodes of east grid for faster plotting
east_line_ylist
list of line nodes of east grid for faster plotting
ew_limits
(min, max) limits of e-w in map_scale units default is None and scales to station area
fig
matplotlib.figure instance for figure
fig_aspect
aspect ratio of plots. default is 1
fig_dpi
resolution of figure in dots-per-inch default is 300
fig_num
figure instance number
fig_size
[width, height] of figure window. default is [6,6]
font_dict
dictionary of font keywords, internally created
font_size
size of ticklables in points, axes labes are font_size+2. default is 4
grid_east
relative location of grid nodes in e-w direction in map_scale units
grid_north
relative location of grid nodes in n-s direction in map_scale units
grid_z
relative location of grid nodes in z direction in map_scale units
index_east
index value of grid_east being plotted
index_north
index value of grid_north being plotted
index_vertical
index value of grid_z being plotted
initial_fn
full path to initial file
key_press
matplotlib.canvas.connect instance
map_scale
[ ‘m’ | ‘km’ ] scale of map. default is km
mesh_east
np.meshgrid(grid_east, grid_north)[0]
mesh_en_east
np.meshgrid(grid_east, grid_north)[0]
mesh_en_north
np.meshgrid(grid_east, grid_north)[1]
mesh_ez_east
np.meshgrid(grid_east, grid_z)[0]
mesh_ez_vertical
np.meshgrid(grid_east, grid_z)[1]
mesh_north
np.meshgrid(grid_east, grid_north)[1]
mesh_nz_north
np.meshgrid(grid_north, grid_z)[0]
mesh_nz_vertical
np.meshgrid(grid_north, grid_z)[1]
model_fn
full path to model file
ms
size of station markers in points. default is 2
nodes_east
relative distance betwen nodes in e-w direction in map_scale units
nodes_north
relative distance betwen nodes in n-s direction in map_scale units
nodes_z
relative distance betwen nodes in z direction in map_scale units
north_line_xlist
list of line nodes north grid for faster plotting
north_line_ylist
list of line nodes north grid for faster plotting
ns_limits
(min, max) limits of plots in n-s direction default is None, set veiwing area to station area
plot_yn
[ ‘y’ | ‘n’ ] ‘y’ to plot on instantiation default is ‘y’
plot_stations
default False
plot_grid
show grid on exported plot; default False
res_model
np.ndarray(n_north, n_east, n_vertical) of model resistivity values in linear scale
save_format
exported format; default png
save_path
path to save exported plots to; default current working folder
station_color
color of station marker. default is black
station_dict_east
location of stations for each east grid row
station_dict_north
location of stations for each north grid row
station_east
location of stations in east direction
station_fn
full path to station file
station_font_color
color of station label
station_font_pad
padding between station marker and label
station_font_rotation
angle of station label
station_font_size
font size of station label
station_font_weight
weight of font for station label
station_id
[min, max] index values for station labels
station_marker
station marker
station_names
name of stations
station_north
location of stations in north direction
subplot_bottom
distance between axes and bottom of figure window
subplot_hspace
distance between subplots in vertical direction
subplot_left
distance between axes and left of figure window
subplot_right
distance between axes and right of figure window
subplot_top
distance between axes and top of figure window
subplot_wspace
distance between subplots in horizontal direction
title
title of plot
xminorticks
location of xminorticks
yminorticks
location of yminorticks
z_limits
(min, max) limits in vertical direction,
Methods
basemap_plot
(depth[, basemap, ...])plot model depth slice on a basemap using basemap modules in matplotlib
export_slices
([plane, indexlist, ...])Plot Slices
get_slice
([option, coords, nsteps, nn, p, ...])- param option
can be either of 'STA', 'XY' or 'XYZ'. For 'STA' or 'XY', a vertical
get the grid line on which a station resides for plotting
on_key_press
(event)on a key press change the slices
plot
()plot:
plot_resistivity_on_seismic
(segy_fn[, ...])- param segy_fn
SegY file name
read in the files to get appropriate information
redraw plot if parameters were changed
save_figure
([save_fn, fig_dpi, file_format, ...])save_figure will save the figure to save_fn.
- basemap_plot(depth, basemap=None, tick_interval=None, save=False, save_path=None, new_figure=True, mesh_rotation_angle=0.0, overlay=False, clip=[0, 0], **basemap_kwargs)[source]¶
plot model depth slice on a basemap using basemap modules in matplotlib
- Parameters
depth – depth in model to plot
tick_interval – tick interval on map in degrees, if None it is calculated from the data extent
save – True/False, whether or not to save and close figure
savepath – full path of file to save to, if None, saves to self.save_path
mesh_rotation_angle – rotation angle of mesh, in degrees clockwise from north
**basemap_kwargs –
provide any valid arguments to Basemap instance and these will be passed to the map.
- New_figure
True/False, whether or not to initiate a new figure for the plot
- export_slices(plane='N-E', indexlist=[], station_buffer=200, save=True)[source]¶
Plot Slices
- Parameters
plane – must be either ‘N-E’, ‘N-Z’ or ‘E-Z’
indexlist – must be a list or 1d numpy array of indices
station_buffer – spatial buffer (in metres) used around grid locations for selecting stations to be projected and plotted on profiles. Ignored if .plot_stations is set to False.
- Returns
[figlist, savepaths]. A list containing (1) lists of Figure objects, for further manipulation (2) corresponding paths for saving them to disk
- get_slice(option='STA', coords=[], nsteps=-1, nn=1, p=4, absolute_query_locations=False, extrapolate=True, reorder_coordinates=False)[source]¶
- Parameters
option – can be either of ‘STA’, ‘XY’ or ‘XYZ’. For ‘STA’ or ‘XY’, a vertical profile is returned based on station coordinates or arbitrary XY coordinates, respectively. For ‘XYZ’, interpolated values at those coordinates are returned
coords – Numpy array of shape (np, 2) for option=’XY’ or of shape (np, 3) for option=’XYZ’, where np is the number of coordinates. Not used for option=’STA’, in which case a vertical profile is created based on station locations.
nsteps – When option is set to ‘STA’ or ‘XY’, nsteps can be used to create a regular grid along the profile in the horizontal direction. By default, when nsteps=-1, the horizontal grid points are defined by station locations or values in the XY array. This parameter is ignored for option=’XYZ’
nn – Number of neighbours to use for interpolation. Nearest neighbour interpolation is returned when nn=1 (default). When nn>1, inverse distance weighted interpolation is returned. See link below for more details: https://en.wikipedia.org/wiki/Inverse_distance_weighting
p – Power parameter, which determines the relative influence of near and far neighbours during interpolation. For p<=3, causes interpolated values to be dominated by points far away. Larger values of p assign greater influence to values near the interpolated point.
absolute_query_locations – if True, query locations are shifted to be centered on the center of station locations; default False, in which case the function treats query locations as relative coordinates. For option=’STA’, this parameter is ignored, since station locations are internally treated as relative coordinates
extrapolate – Extrapolates values (default), which can be particularly useful for extracting values at nodes, since the field values are given for cell-centres.
reorder_coordinates – attempts to reorder coordinates (when option is ‘STA’ or ‘XY’) to form a continuous line.
- Returns
- 1: when option is ‘STA’ or ‘XY’
gd, gz, gv : where gd, gz and gv are 2D grids of distance (along profile), depth and interpolated values, respectively. The shape of the 2D grids depend on the number of stations or the number of xy coordinates provided, for options ‘STA’ or ‘XY’, respectively, the number of vertical model grid points and whether regular gridding in the horizontal direction was enabled with nsteps>-1.
- 2: when option is ‘XYZ’
gv : list of interpolated values of shape (np)
- plot_resistivity_on_seismic(segy_fn, velocity_model=6000, pick_every=10, ax=None, cb_ax=None, percent_clip=99, alpha=0.5, **kwargs)[source]¶
- Parameters
segy_fn – SegY file name
velocity_model – can be either the name of a velocity-model file containing stacking velocities for the given 2D seismic line, or a floating point value representing a constant velocity (m/s)
pick_every – this parameter controls the decimation factor for the SegY file; e.g. if pick_every=10, every 10th trace from the SegY file is read in. This significantly speeds up plotting routines.
ax – figure axes
cb_ax – colorbar axes
percent_clip – percentile value used for filtering out seismic amplitudes from plot; e.g. for a value of 99, only seismic amplitudes above the 99th percentile are plotted. The parameter is tuned to plot only the required level of seismic detail.
alpha – alpha value used while resistivity and seismic values
kwargs –
max_depth : maximum depth extent of plots time_shift : time shift in ms to remove topography
- Returns
fig, ax : a figure and an plot axes object are returned when the parameter ax is not provided
- 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.plotAllResponses() >>> #change line width >>> p1.lw = 2 >>> p1.redraw_plot()
- class mtpy.modeling.modem.Residual(**kwargs)[source]¶
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
Methods
calculate_residual_from_data
([data_fn, ...])created by ak on 26/09/2017
write rms station data to file
get_rms
read_residual_file
- class mtpy.modeling.modem.Stations(**kwargs)[source]¶
station locations class
- ..note:: If the survey steps across multiple UTM zones, then a
distance will be added to the stations to place them in the correct location. This distance is _utm_grid_size_north and _utm_grid_size_east. You should these parameters to place the locations in the proper spot as grid distances and overlaps change over the globe. This is not implemented yet
- Attributes
center_point
calculate the center point from the given station locations
- east
- elev
- lat
- lon
- north
- rel_east
- rel_elev
- rel_north
- station
- utm_zone
Methods
calculate_rel_locations
([shift_east, ...])put station in a coordinate system relative to (shift_east, shift_north) (+) shift right or up (-) shift left or down
If the stations cross utm zones, then estimate distance by computing distance on a sphere.
get_station_locations
(input_list)get station locations from a list of edi files
rotate_stations
(rotation_angle)Rotate stations assuming N is 0
- 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
- Returns
- **center_location**np.ndarray
structured array of length 1 dtype includes (east, north, zone, lat, lon)
- check_utm_crossing()[source]¶
If the stations cross utm zones, then estimate distance by computing distance on a sphere.
# Generate files for ModEM
# revised by JP 2017 # revised by AK 2017 to bring across functionality from ak branch
- class mtpy.modeling.modem.plot_response.PlotResponse(data_fn=None, resp_fn=None, **kwargs)[source]¶
plot data and response
Plots the real and imaginary impedance and induction vector if present.
- Example
>>> import mtpy.modeling.modem as modem >>> dfn = r"/home/MT/ModEM/Inv1/DataFile.dat" >>> rfn = r"/home/MT/ModEM/Inv1/Test_resp_000.dat" >>> mrp = modem.PlotResponse(data_fn=dfn, resp_fn=rfn) >>> # plot only the TE and TM modes >>> mrp.plot_component = 2 >>> mrp.redraw_plot()
Attributes
Description
color_mode
[ ‘color’ | ‘bw’ ] color or black and white plots
cted
color for data Z_XX and Z_XY mode
ctem
color for model Z_XX and Z_XY mode
ctmd
color for data Z_YX and Z_YY mode
ctmm
color for model Z_YX and Z_YY mode
data_fn
full path to data file
data_object
WSResponse instance
e_capsize
cap size of error bars in points (default is .5)
e_capthick
cap thickness of error bars in points (default is 1)
fig_dpi
resolution of figure in dots-per-inch (300)
fig_list
list of matplotlib.figure instances for plots
fig_size
size of figure in inches (default is [6, 6])
font_size
size of font for tick labels, axes labels are font_size+2 (default is 7)
legend_border_axes_pad
padding between legend box and axes
legend_border_pad
padding between border of legend and symbols
legend_handle_text_pad
padding between text labels and symbols of legend
legend_label_spacing
padding between labels
legend_loc
location of legend
legend_marker_scale
scale of symbols in legend
lw
line width data curves (default is .5)
ms
size of markers (default is 1.5)
lw_r
line width response curves (default is .5)
ms_r
size of markers response curves (default is 1.5)
mted
marker for data Z_XX and Z_XY mode
mtem
marker for model Z_XX and Z_XY mode
mtmd
marker for data Z_YX and Z_YY mode
mtmm
marker for model Z_YX and Z_YY mode
phase_limits
limits of phase
plot_component
[ 2 | 4 ] 2 for TE and TM or 4 for all components
plot_style
[ 1 | 2 ] 1 to plot each mode in a seperate subplot and 2 to plot xx, xy and yx, yy in same plots
plot_type
[ ‘1’ | list of station name ] ‘1’ to plot all stations in data file or input a list of station names to plot if station_fn is input, otherwise input a list of integers associated with the index with in the data file, ie 2 for 2nd station
plot_z
[ True | False ] default is True to plot impedance, False for plotting resistivity and phase
plot_yn
[ ‘n’ | ‘y’ ] to plot on instantiation
res_limits
limits of resistivity in linear scale
resp_fn
full path to response file
resp_object
WSResponse object for resp_fn, or list of WSResponse objects if resp_fn is a list of response files
station_fn
full path to station file written by WSStation
subplot_bottom
space between axes and bottom of figure
subplot_hspace
space between subplots in vertical direction
subplot_left
space between axes and left of figure
subplot_right
space between axes and right of figure
subplot_top
space between axes and top of figure
subplot_wspace
space between subplots in horizontal direction
Methods
redraw plot if parameters were changed
save_figure
(save_fn[, file_format, ...])save_plot will save the figure to save_fn.
plot
- 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.plotAllResponses() >>> #change line width >>> p1.lw = 2 >>> p1.redraw_plot()
# Generate files for ModEM
# revised by JP 2017 # revised by AK 2017 to bring across functionality from ak branch
- class mtpy.modeling.modem.plot_slices.PlotSlices(model_fn, data_fn=None, **kwargs)[source]¶
Plot all cartesian axis-aligned slices and be able to scroll through the model
Extract arbitrary profiles (e.g. along a seismic line) from a model
- Example
>>> import mtpy.modeling.modem as modem >>> mfn = r"/home/modem/Inv1/Modular_NLCG_100.rho" >>> dfn = r"/home/modem/Inv1/ModEM_data.dat" >>> pds = ws.PlotSlices(model_fn=mfn, data_fn=dfn)
Buttons
Description
‘e’
moves n-s slice east by one model block
‘w’
moves n-s slice west by one model block
‘n’
moves e-w slice north by one model block
‘m’
moves e-w slice south by one model block
‘d’
moves depth slice down by one model block
‘u’
moves depth slice up by one model block
Attributes
Description
ax_en
matplotlib.axes instance for depth slice map view
ax_ez
matplotlib.axes instance for e-w slice
ax_map
matplotlib.axes instance for location map
ax_nz
matplotlib.axes instance for n-s slice
climits
(min , max) color limits on resistivity in log scale. default is (0, 4)
cmap
name of color map for resisitiviy. default is ‘jet_r’
data_fn
full path to data file name
draw_colorbar
show colorbar on exported plot; default True
dscale
scaling parameter depending on map_scale
east_line_xlist
list of line nodes of east grid for faster plotting
east_line_ylist
list of line nodes of east grid for faster plotting
ew_limits
(min, max) limits of e-w in map_scale units default is None and scales to station area
fig
matplotlib.figure instance for figure
fig_aspect
aspect ratio of plots. default is 1
fig_dpi
resolution of figure in dots-per-inch default is 300
fig_num
figure instance number
fig_size
[width, height] of figure window. default is [6,6]
font_dict
dictionary of font keywords, internally created
font_size
size of ticklables in points, axes labes are font_size+2. default is 4
grid_east
relative location of grid nodes in e-w direction in map_scale units
grid_north
relative location of grid nodes in n-s direction in map_scale units
grid_z
relative location of grid nodes in z direction in map_scale units
index_east
index value of grid_east being plotted
index_north
index value of grid_north being plotted
index_vertical
index value of grid_z being plotted
initial_fn
full path to initial file
key_press
matplotlib.canvas.connect instance
map_scale
[ ‘m’ | ‘km’ ] scale of map. default is km
mesh_east
np.meshgrid(grid_east, grid_north)[0]
mesh_en_east
np.meshgrid(grid_east, grid_north)[0]
mesh_en_north
np.meshgrid(grid_east, grid_north)[1]
mesh_ez_east
np.meshgrid(grid_east, grid_z)[0]
mesh_ez_vertical
np.meshgrid(grid_east, grid_z)[1]
mesh_north
np.meshgrid(grid_east, grid_north)[1]
mesh_nz_north
np.meshgrid(grid_north, grid_z)[0]
mesh_nz_vertical
np.meshgrid(grid_north, grid_z)[1]
model_fn
full path to model file
ms
size of station markers in points. default is 2
nodes_east
relative distance betwen nodes in e-w direction in map_scale units
nodes_north
relative distance betwen nodes in n-s direction in map_scale units
nodes_z
relative distance betwen nodes in z direction in map_scale units
north_line_xlist
list of line nodes north grid for faster plotting
north_line_ylist
list of line nodes north grid for faster plotting
ns_limits
(min, max) limits of plots in n-s direction default is None, set veiwing area to station area
plot_yn
[ ‘y’ | ‘n’ ] ‘y’ to plot on instantiation default is ‘y’
plot_stations
default False
plot_grid
show grid on exported plot; default False
res_model
np.ndarray(n_north, n_east, n_vertical) of model resistivity values in linear scale
save_format
exported format; default png
save_path
path to save exported plots to; default current working folder
station_color
color of station marker. default is black
station_dict_east
location of stations for each east grid row
station_dict_north
location of stations for each north grid row
station_east
location of stations in east direction
station_fn
full path to station file
station_font_color
color of station label
station_font_pad
padding between station marker and label
station_font_rotation
angle of station label
station_font_size
font size of station label
station_font_weight
weight of font for station label
station_id
[min, max] index values for station labels
station_marker
station marker
station_names
name of stations
station_north
location of stations in north direction
subplot_bottom
distance between axes and bottom of figure window
subplot_hspace
distance between subplots in vertical direction
subplot_left
distance between axes and left of figure window
subplot_right
distance between axes and right of figure window
subplot_top
distance between axes and top of figure window
subplot_wspace
distance between subplots in horizontal direction
title
title of plot
xminorticks
location of xminorticks
yminorticks
location of yminorticks
z_limits
(min, max) limits in vertical direction,
Methods
basemap_plot
(depth[, basemap, ...])plot model depth slice on a basemap using basemap modules in matplotlib
export_slices
([plane, indexlist, ...])Plot Slices
get_slice
([option, coords, nsteps, nn, p, ...])- param option
can be either of 'STA', 'XY' or 'XYZ'. For 'STA' or 'XY', a vertical
get the grid line on which a station resides for plotting
on_key_press
(event)on a key press change the slices
plot
()plot:
plot_resistivity_on_seismic
(segy_fn[, ...])- param segy_fn
SegY file name
read in the files to get appropriate information
redraw plot if parameters were changed
save_figure
([save_fn, fig_dpi, file_format, ...])save_figure will save the figure to save_fn.
- basemap_plot(depth, basemap=None, tick_interval=None, save=False, save_path=None, new_figure=True, mesh_rotation_angle=0.0, overlay=False, clip=[0, 0], **basemap_kwargs)[source]¶
plot model depth slice on a basemap using basemap modules in matplotlib
- Parameters
depth – depth in model to plot
tick_interval – tick interval on map in degrees, if None it is calculated from the data extent
save – True/False, whether or not to save and close figure
savepath – full path of file to save to, if None, saves to self.save_path
mesh_rotation_angle – rotation angle of mesh, in degrees clockwise from north
**basemap_kwargs –
provide any valid arguments to Basemap instance and these will be passed to the map.
- New_figure
True/False, whether or not to initiate a new figure for the plot
- export_slices(plane='N-E', indexlist=[], station_buffer=200, save=True)[source]¶
Plot Slices
- Parameters
plane – must be either ‘N-E’, ‘N-Z’ or ‘E-Z’
indexlist – must be a list or 1d numpy array of indices
station_buffer – spatial buffer (in metres) used around grid locations for selecting stations to be projected and plotted on profiles. Ignored if .plot_stations is set to False.
- Returns
[figlist, savepaths]. A list containing (1) lists of Figure objects, for further manipulation (2) corresponding paths for saving them to disk
- get_slice(option='STA', coords=[], nsteps=-1, nn=1, p=4, absolute_query_locations=False, extrapolate=True, reorder_coordinates=False)[source]¶
- Parameters
option – can be either of ‘STA’, ‘XY’ or ‘XYZ’. For ‘STA’ or ‘XY’, a vertical profile is returned based on station coordinates or arbitrary XY coordinates, respectively. For ‘XYZ’, interpolated values at those coordinates are returned
coords – Numpy array of shape (np, 2) for option=’XY’ or of shape (np, 3) for option=’XYZ’, where np is the number of coordinates. Not used for option=’STA’, in which case a vertical profile is created based on station locations.
nsteps – When option is set to ‘STA’ or ‘XY’, nsteps can be used to create a regular grid along the profile in the horizontal direction. By default, when nsteps=-1, the horizontal grid points are defined by station locations or values in the XY array. This parameter is ignored for option=’XYZ’
nn – Number of neighbours to use for interpolation. Nearest neighbour interpolation is returned when nn=1 (default). When nn>1, inverse distance weighted interpolation is returned. See link below for more details: https://en.wikipedia.org/wiki/Inverse_distance_weighting
p – Power parameter, which determines the relative influence of near and far neighbours during interpolation. For p<=3, causes interpolated values to be dominated by points far away. Larger values of p assign greater influence to values near the interpolated point.
absolute_query_locations – if True, query locations are shifted to be centered on the center of station locations; default False, in which case the function treats query locations as relative coordinates. For option=’STA’, this parameter is ignored, since station locations are internally treated as relative coordinates
extrapolate – Extrapolates values (default), which can be particularly useful for extracting values at nodes, since the field values are given for cell-centres.
reorder_coordinates – attempts to reorder coordinates (when option is ‘STA’ or ‘XY’) to form a continuous line.
- Returns
- 1: when option is ‘STA’ or ‘XY’
gd, gz, gv : where gd, gz and gv are 2D grids of distance (along profile), depth and interpolated values, respectively. The shape of the 2D grids depend on the number of stations or the number of xy coordinates provided, for options ‘STA’ or ‘XY’, respectively, the number of vertical model grid points and whether regular gridding in the horizontal direction was enabled with nsteps>-1.
- 2: when option is ‘XYZ’
gv : list of interpolated values of shape (np)
- plot_resistivity_on_seismic(segy_fn, velocity_model=6000, pick_every=10, ax=None, cb_ax=None, percent_clip=99, alpha=0.5, **kwargs)[source]¶
- Parameters
segy_fn – SegY file name
velocity_model – can be either the name of a velocity-model file containing stacking velocities for the given 2D seismic line, or a floating point value representing a constant velocity (m/s)
pick_every – this parameter controls the decimation factor for the SegY file; e.g. if pick_every=10, every 10th trace from the SegY file is read in. This significantly speeds up plotting routines.
ax – figure axes
cb_ax – colorbar axes
percent_clip – percentile value used for filtering out seismic amplitudes from plot; e.g. for a value of 99, only seismic amplitudes above the 99th percentile are plotted. The parameter is tuned to plot only the required level of seismic detail.
alpha – alpha value used while resistivity and seismic values
kwargs –
max_depth : maximum depth extent of plots time_shift : time shift in ms to remove topography
- Returns
fig, ax : a figure and an plot axes object are returned when the parameter ax is not provided
- 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.plotAllResponses() >>> #change line width >>> p1.lw = 2 >>> p1.redraw_plot()
Create Phase Tensor Map from the ModEM’s output Resistivity model
- class mtpy.modeling.modem.phase_tensor_maps.PlotPTMaps(data_fn=None, resp_fn=None, model_fn=None, **kwargs)[source]¶
Plot phase tensor maps including residual pt if response file is input.
- Plot only data for one period
>>> import mtpy.modeling.ws3dinv as ws >>> dfn = r"/home/MT/ws3dinv/Inv1/WSDataFile.dat" >>> ptm = ws.PlotPTMaps(data_fn=dfn, plot_period_list=[0])
- Plot data and model response
>>> import mtpy.modeling.ws3dinv as ws >>> dfn = r"/home/MT/ws3dinv/Inv1/WSDataFile.dat" >>> rfn = r"/home/MT/ws3dinv/Inv1/Test_resp.00" >>> mfn = r"/home/MT/ws3dinv/Inv1/Test_model.00" >>> ptm = ws.PlotPTMaps(data_fn=dfn, resp_fn=rfn, model_fn=mfn, >>> ... plot_period_list=[0]) >>> # adjust colorbar >>> ptm.cb_res_pad = 1.25 >>> ptm.redraw_plot()
Methods
get_period_attributes
(periodIdx, key[, ptarray])Returns, for a given period, a list of attribute values for key (e.g.
plot
([period, periodIdx, save2file])Plot phase tensor maps for data and or response, each figure is of a different period.
plot_on_axes
(ax, m, periodIdx[, ptarray, ...])Plots phase tensors for a given period index.
redraw plot if parameters were changed
save_all_figures
([save_path, fig_dpi, ...])save_figure will save all figures in fig_list to save_fn.
write_pt_data_to_gmt
([period, epsg, ...])write data to plot phase tensor ellipses in gmt.
write_pt_data_to_text
- get_period_attributes(periodIdx, key, ptarray='data')[source]¶
Returns, for a given period, a list of attribute values for key (e.g. skew, phimax, etc.).
- Parameters
periodIdx – index of period; print out _plot_period for periods available
key – attribute key
ptarray – name of data-array to access for retrieving attributes; can be either ‘data’, ‘resp’ or ‘resid’
- Returns
numpy array of attribute values
- plot(period=None, periodIdx=0, save2file=None, **kwargs)[source]¶
Plot phase tensor maps for data and or response, each figure is of a different period. If response is input a third column is added which is the residual phase tensor showing where the model is not fitting the data well. The data is plotted in km.
- Args:
period: the period index to plot, default=0
Returns:
- plot_on_axes(ax, m, periodIdx, ptarray='data', ellipse_size_factor=10000, cvals=None, map_scale='m', centre_shift=[0, 0], plot_tipper='n', tipper_size_factor=100000.0, **kwargs)[source]¶
Plots phase tensors for a given period index.
- Parameters
ax – plot axis
m – basemap instance
periodIdx – period index
ptarray – name of data-array to access for retrieving attributes; can be either ‘data’, ‘resp’ or ‘resid’
ellipse_size_factor – factor to control ellipse size
cvals – list of colour values for colouring each ellipse; must be of the same length as the number of tuples for each period
map_scale – map length scale
kwargs – list of relevant matplotlib arguments (e.g. zorder, alpha, etc.)
plot_tipper – string (‘n’, ‘yr’, ‘yi’, or ‘yri’) to plot no tipper, real only, imaginary only, or both
tipper_size_factor – scaling factor for tipper vectors
- 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.plotAllResponses() >>> #change line width >>> p1.lw = 2 >>> p1.redraw_plot()
- save_all_figures(save_path=None, fig_dpi=None, file_format='pdf', orientation='landscape', close_fig='y')[source]¶
save_figure will save all figures in fig_list to save_fn.
- write_pt_data_to_gmt(period=None, epsg=None, savepath='.', center_utm=None, colorby='phimin', attribute='data', clim=None)[source]¶
write data to plot phase tensor ellipses in gmt. saves a gmt script and text file containing ellipse data
provide: period to plot (seconds) epsg for the projection the model was projected to (google “epsg your_projection_name” and you will find it) centre_utm - utm coordinates for centre position of model, if not
provided, script will try and extract it from data file
colorby - what to colour the ellipses by, ‘phimin’, ‘phimax’, or ‘skew’ attribute - attribute to plot ‘data’, ‘resp’, or ‘resid’ for data,
response or residuals
Module Occam 1D¶
- Wrapper class to interact with Occam1D written by Kerry Keys at Scripps
adapted from the method of Constable et al., [1987].
This class only deals with the MT functionality of the Fortran code, so it can make the input files for computing the 1D MT response of an input model and or data. It can also read the output and plot them in a useful way.
Note that when you run the inversion code, the convergence is quite quick, within the first few iterations, so have a look at the L2 cure to decide which iteration to plot, otherwise if you look at iterations long after convergence the models will be unreliable.
Key, K., 2009, 1D inversion of multicomponent, multi-frequency marine CSEM data: Methodology and synthetic studies for resolving thin resistive layers: Geophysics, 74, F9–F20.
The original paper describing the Occam’s inversion approach is:
Constable, S. C., R. L. Parker, and C. G. Constable, 1987, Occam’s inversion –– A practical algorithm for generating smooth models from electromagnetic sounding data, Geophysics, 52 (03), 289–300.
- Intended Use
>>> import mtpy.modeling.occam1d as occam1d >>> #--> make a data file >>> d1 = occam1d.Data() >>> d1.write_data_file(edi_file=r'/home/MT/mt01.edi', res_err=10, phase_err=2.5, >>> ... save_path=r"/home/occam1d/mt01/TE", mode='TE') >>> #--> make a model file >>> m1 = occam1d.Model() >>> m1.write_model_file(save_path=d1.save_path, target_depth=15000) >>> #--> make a startup file >>> s1 = occam1d.Startup() >>> s1.data_fn = d1.data_fn >>> s1.model_fn = m1.model_fn >>> s1.save_path = m1.save_path >>> s1.write_startup_file() >>> #--> run occam1d from python >>> occam_path = r"/home/occam1d/Occam1D_executable" >>> occam1d.Run(s1.startup_fn, occam_path, mode='TE') >>> #--plot the L2 curve >>> l2 = occam1d.PlotL2(d1.save_path, m1.model_fn) >>> #--> see that iteration 7 is the optimum model to plot >>> p1 = occam1d.Plot1DResponse() >>> p1.data_te_fn = d1.data_fn >>> p1.model_fn = m1.model_fn >>> p1.iter_te_fn = r"/home/occam1d/mt01/TE/TE_7.iter" >>> p1.resp_te_fn = r"/home/occam1d/mt01/TE/TE_7.resp" >>> p1.plot()
@author: J. Peacock (Oct. 2013)
- class mtpy.modeling.occam1d.Data(data_fn=None, **kwargs)[source]¶
reads and writes occam 1D data files
Attributes
Description
_data_fn
basename of data file default is Occam1DDataFile
_header_line
header line for description of data columns
_ss
string spacing default is 6*’ ‘
_string_fmt
format of data default is ‘+.6e’
data
array of data
data_fn
full path to data file
freq
frequency array of data
mode
mode to invert for [ ‘TE’ | ‘TM’ | ‘det’ ]
phase_te
array of TE phase
phase_tm
array of TM phase
res_te
array of TE apparent resistivity
res_tm
array of TM apparent resistivity
resp_fn
full path to response file
save_path
path to save files to
Methods
Description
write_data_file
write an Occam1D data file
read_data_file
read an Occam1D data file
read_resp_file
read a .resp file output by Occam1D
- Example
>>> import mtpy.modeling.occam1d as occam1d >>> #--> make a data file for TE mode >>> d1 = occam1d.Data() >>> d1.write_data_file(edi_file=r'/home/MT/mt01.edi', res_err=10, phase_err=2.5, >>> ... save_path=r"/home/occam1d/mt01/TE", mode='TE')
Methods
read_data_file
([data_fn])reads a 1D data file
read_resp_file
([resp_fn, data_fn])read response file
write_data_file
([rp_tuple, edi_file, ...])make1Ddatafile will write a data file for Occam1D
- read_resp_file(resp_fn=None, data_fn=None)[source]¶
read response file
resp_fn : full path to response file
data_fn : full path to data file
- write_data_file(rp_tuple=None, edi_file=None, save_path=None, mode='det', res_err='data', phase_err='data', thetar=0, res_errorfloor=0.0, phase_errorfloor=0.0, z_errorfloor=0.0, remove_outofquadrant=False)[source]¶
make1Ddatafile will write a data file for Occam1D
Arguments:¶
- rp_tuplenp.ndarray (freq, res, res_err, phase, phase_err)
with res, phase having shape (num_freq, 2, 2).
- edi_filestring
full path to edi file to be modeled.
- save_pathstring
path to save the file, if None set to dirname of station if edipath = None. Otherwise set to dirname of edipath.
- thetarfloat
rotation angle to rotate Z. Clockwise positive and N=0 default = 0
- mode[ ‘te’ | ‘tm’ | ‘det’]
- mode to model can be (*default*=’both’):
‘te’ for just TE mode (res/phase)
‘tm’ for just TM mode (res/phase)
- ‘det’ for the determinant of Z (converted to
res/phase)
add ‘z’ to any of these options to model impedance tensor values instead of res/phase
- res_errfloat
errorbar for resistivity values. Can be set to ( default = ‘data’):
‘data’ for errorbars from the data
percent number ex. 10 for ten percent
- phase_errfloat
errorbar for phase values. Can be set to ( default = ‘data’):
‘data’ for errorbars from the data
percent number ex. 10 for ten percent
- res_errorfloor: float
error floor for resistivity values in percent
- phase_errorfloor: float
error floor for phase in degrees
- remove_outofquadrant: True/False; option to remove the resistivity and
phase values for points with phases out of the 1st/3rd quadrant (occam requires 0 < phase < 90 degrees; phases in the 3rd quadrant are shifted to the first by adding 180 degrees)
- Example
>>> import mtpy.modeling.occam1d as occam1d >>> #--> make a data file >>> d1 = occam1d.Data() >>> d1.write_data_file(edi_file=r'/home/MT/mt01.edi', res_err=10, >>> ... phase_err=2.5, mode='TE', >>> ... save_path=r"/home/occam1d/mt01/TE")
- class mtpy.modeling.occam1d.Model(model_fn=None, **kwargs)[source]¶
read and write the model file fo Occam1D
All depth measurements are in meters.
Attributes
Description
_model_fn
basename for model file default is Model1D
_ss
string spacing in model file default is 3*’ ‘
_string_fmt
format of model layers default is ‘.0f’
air_layer_height
height of air layer default is 10000
bottom_layer
bottom of the model default is 50000
itdict
dictionary of values from iteration file
iter_fn
full path to iteration file
model_depth
array of model depths
model_fn
full path to model file
model_penalty
array of penalties for each model layer
model_preference_penalty
array of model preference penalties for each layer
model_prefernce
array of preferences for each layer
model_res
array of resistivities for each layer
n_layers
number of layers in the model
num_params
number of parameters to invert for (n_layers+2)
pad_z
padding of model at depth default is 5 blocks
save_path
path to save files
target_depth
depth of target to investigate
z1_layer
depth of first layer default is 10
Methods
Description
write_model_file
write an Occam1D model file, where depth increases on a logarithmic scale
read_model_file
read an Occam1D model file
read_iter_file
read an .iter file output by Occam1D
- Example
>>> #--> make a model file >>> m1 = occam1d.Model() >>> m1.write_model_file(save_path=r"/home/occam1d/mt01/TE")
Methods
read_iter_file
([iter_fn, model_fn])read an 1D iteration file
read_model_file
([model_fn])will read in model 1D file
write_model_file
([save_path])Makes a 1D model file for Occam1D.
- class mtpy.modeling.occam1d.Plot1DResponse(data_te_fn=None, data_tm_fn=None, model_fn=None, resp_te_fn=None, resp_tm_fn=None, iter_te_fn=None, iter_tm_fn=None, **kwargs)[source]¶
plot the 1D response and model. Plots apparent resisitivity and phase in different subplots with the model on the far right. You can plot both TE and TM modes together along with different iterations of the model. These will be plotted in different colors or shades of gray depneng on color_scale.
- Example
>>> import mtpy.modeling.occam1d as occam1d >>> p1 = occam1d.Plot1DResponse(plot_yn='n') >>> p1.data_te_fn = r"/home/occam1d/mt01/TE/Occam_DataFile_TE.dat" >>> p1.data_tm_fn = r"/home/occam1d/mt01/TM/Occam_DataFile_TM.dat" >>> p1.model_fn = r"/home/occam1d/mt01/TE/Model1D" >>> p1.iter_te_fn = [r"/home/occam1d/mt01/TE/TE_{0}.iter".format(ii) >>> ... for ii in range(5,10)] >>> p1.iter_tm_fn = [r"/home/occam1d/mt01/TM/TM_{0}.iter".format(ii) >>> ... for ii in range(5,10)] >>> p1.resp_te_fn = [r"/home/occam1d/mt01/TE/TE_{0}.resp".format(ii) >>> ... for ii in range(5,10)] >>> p1.resp_tm_fn = [r"/home/occam1d/mt01/TM/TM_{0}.resp".format(ii) >>> ... for ii in range(5,10)] >>> p1.plot()
Attributes
Description
axm
matplotlib.axes instance for model subplot
axp
matplotlib.axes instance for phase subplot
axr
matplotlib.axes instance for app. res subplot
color_mode
[ ‘color’ | ‘bw’ ]
cted
color of TE data markers
ctem
color of TM data markers
ctmd
color of TE model markers
ctmm
color of TM model markers
data_te_fn
full path to data file for TE mode
data_tm_fn
full path to data file for TM mode
depth_limits
(min, max) limits for depth plot in depth_units
depth_scale
[ ‘log’ | ‘linear’ ] default is linear
depth_units
[ ‘m’ | ‘km’ ] *default is ‘km’
e_capsize
capsize of error bars
e_capthick
cap thickness of error bars
fig
matplotlib.figure instance for plot
fig_dpi
resolution in dots-per-inch for figure
fig_num
number of figure instance
fig_size
size of figure in inches [width, height]
font_size
size of axes tick labels, axes labels are +2
grid_alpha
transparency of grid
grid_color
color of grid
iter_te_fn
full path or list of .iter files for TE mode
iter_tm_fn
full path or list of .iter files for TM mode
lw
width of lines for model
model_fn
full path to model file
ms
marker size
mted
marker for TE data
mtem
marker for TM data
mtmd
marker for TE model
mtmm
marker for TM model
phase_limits
(min, max) limits on phase in degrees
phase_major_ticks
spacing for major ticks in phase
phase_minor_ticks
spacing for minor ticks in phase
plot_yn
[ ‘y’ | ‘n’ ] plot on instantiation
res_limits
limits of resistivity in linear scale
resp_te_fn
full path or list of .resp files for TE mode
resp_tm_fn
full path or list of .iter files for TM mode
subplot_bottom
spacing of subplots from bottom of figure
subplot_hspace
height spacing between subplots
subplot_left
spacing of subplots from left of figure
subplot_right
spacing of subplots from right of figure
subplot_top
spacing of subplots from top of figure
subplot_wspace
width spacing between subplots
title_str
title of plot
Methods
plot
()plot data, response and model
redraw plot if parameters were changed
save_figure
(save_fn[, file_format, ...])save_plot will save the figure to save_fn.
update_plot
(fig)update any parameters that where changed using the built-in draw from canvas.
- 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.plotAllResponses() >>> #change line width >>> p1.lw = 2 >>> 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.
- update_plot(fig)[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.plotAllResponses() >>> [ax.grid(True, which='major') for ax in [ps1.axrte,ps1.axtep]] >>> ps1.update_plot()
- class mtpy.modeling.occam1d.PlotL2(dir_path, model_fn, **kwargs)[source]¶
plot L2 curve of iteration vs rms and roughness
Methods
plot
()plot L2 curve
redraw plot if parameters were changed
save_figure
(save_fn[, file_format, ...])save_plot will save the figure to save_fn.
update any parameters that where changed using the built-in draw from canvas.
- 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.plotAllResponses() >>> #change line width >>> p1.lw = 2 >>> p1.redraw_plot()
- save_figure(save_fn, file_format='pdf', orientation='portrait', fig_dpi=None, close_fig='y')[source]¶
save_plot will save the figure to save_fn.
- 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.plotAllResponses() >>> [ax.grid(True, which='major') for ax in [ps1.axrte,ps1.axtep]] >>> ps1.update_plot()
- class mtpy.modeling.occam1d.Run(startup_fn=None, occam_path=None, **kwargs)[source]¶
run occam 1d from python given the correct files and location of occam1d executable
Methods
run_occam1d
- class mtpy.modeling.occam1d.Startup(data_fn=None, model_fn=None, **kwargs)[source]¶
read and write input files for Occam1D
Attributes
Description
_ss
string spacing
_startup_fn
basename of startup file default is OccamStartup1D
data_fn
full path to data file
debug_level
debug level default is 1
description
description of inversion for your self default is 1D_Occam_Inv
max_iter
maximum number of iterations default is 20
model_fn
full path to model file
rough_type
roughness type default is 1
save_path
full path to save files to
start_iter
first iteration number default is 0
start_lagrange
starting lagrange number on log scale default is 5
start_misfit
starting misfit value default is 100
start_rho
starting resistivity value (halfspace) in log scale default is 100
start_rough
starting roughness (ignored by Occam1D) default is 1E7
startup_fn
full path to startup file
target_rms
target rms default is 1.0
Methods
read_startup_file
(startup_fn)reads in a 1D input file
write_startup_file
([save_path])Make a 1D input file for Occam 1D
- read_startup_file(startup_fn)[source]¶
reads in a 1D input file
Arguments:¶
inputfn : full path to input file
- write_startup_file(save_path=None, **kwargs)[source]¶
Make a 1D input file for Occam 1D
Arguments:¶
- savepathfull path to save input file to, if just path then
saved as savepath/input
- model_fnfull path to model file, if None then assumed to be in
savepath/model.mod
- data_fnfull path to data file, if None then assumed to be
in savepath/TE.dat or TM.dat
rough_type : roughness type. default = 0
max_iter : maximum number of iterations. default = 20
target_rms : target rms value. default = 1.0
- start_rhostarting resistivity value on linear scale.
default = 100
description : description of the inversion.
- start_lagrangestarting Lagrange multiplier for smoothness.
default = 5
start_rough : starting roughness value. default = 1E7
- debuglevelsomething to do with how Fortran debuggs the code
Almost always leave at default = 1
- start_iterthe starting iteration number, handy if the
starting model is from a previous run. default = 0
start_misfit : starting misfit value. default = 100
- mtpy.modeling.occam1d.build_run()[source]¶
build input files and run a suite of models in series (pretty quick so won’t bother parallelise)
run Occam1d on each set of inputs. Occam is run twice. First to get the lowest possible misfit. we then set the target rms to a factor (default 1.05) times the minimum rms achieved and run to get the smoothest model.
author: Alison Kirkby (2016)
- mtpy.modeling.occam1d.divide_inputs(work_to_do, size)[source]¶
divide list of inputs into chunks to send to each processor
- mtpy.modeling.occam1d.generate_inputfiles(**input_parameters)[source]¶
generate all the input files to run occam1d, return the path and the startup files to run.
author: Alison Kirkby (2016)
- mtpy.modeling.occam1d.get_strike(mt_object, fmin, fmax, strike_approx=0)[source]¶
get the strike from the z array, choosing the strike angle that is closest to the azimuth of the PT ellipse (PT strike).
if there is not strike available from the z array use the PT strike.
Module Occam 2D¶
Spin-off from ‘occamtools’ (Created August 2011, re-written August 2013)
Tools for Occam2D
authors: JP/LK
- Classes:
Data
Model
Setup
Run
Plot
Mask
- Functions:
getdatetime
makestartfiles
writemeshfile
writemodelfile
writestartupfile
read_datafile
get_model_setup
blocks_elements_setup
- class mtpy.modeling.occam2d_rewrite.Data(edi_path=None, **kwargs)[source]¶
Reads and writes data files and more.
Inherets Profile, so the intended use is to use Data to project stations onto a profile, then write the data file.
Model Modes
Description
1 or log_all
Log resistivity of TE and TM plus Tipper
2 or log_te_tip
Log resistivity of TE plus Tipper
3 or log_tm_tip
Log resistivity of TM plus Tipper
4 or log_te_tm
Log resistivity of TE and TM
5 or log_te
Log resistivity of TE
6 or log_tm
Log resistivity of TM
7 or all
TE, TM and Tipper
8 or te_tip
TE plus Tipper
9 or tm_tip
TM plus Tipper
10 or te_tm
TE and TM mode
11 or te
TE mode
12 or tm
TM mode
13 or tip
Only Tipper
- datais a list of dictioinaries containing the data for each station.
- keys include:
‘station’ – name of station
‘offset’ – profile line offset
‘te_res’ – TE resisitivity in linear scale
‘tm_res’ – TM resistivity in linear scale
‘te_phase’ – TE phase in degrees
‘tm_phase’ – TM phase in degrees in first quadrant
‘re_tip’ – real part of tipper along profile
‘im_tip’ – imaginary part of tipper along profile
each key is a np.ndarray(2, num_freq) index 0 is for data index 1 is for error
Key Words/Attributes
Desctription
_data_header
header line in data file
_data_string
full data string
_profile_generated
[ True | False ] True if profile has already been generated.
_rotate_to_strike
[ True | False ] True to rotate data to strike angle. default is True
data
list of dictionaries of data for each station. see above
data_fn
full path to data file
data_list
list of lines to write to data file
edi_list
list of mtpy.core.mt instances for each .edi file read
edi_path
directory path where .edi files are
edi_type
[ ‘z’ | ‘spectra’ ] for .edi format
elevation_model
model elevation np.ndarray(east, north, elevation) in meters
elevation_profile
elevation along profile np.ndarray (x, elev) (m)
fn_basename
data file basename default is OccamDataFile.dat
freq
list of frequencies to use for the inversion
freq_max
max frequency to use in inversion. default is None
freq_min
min frequency to use in inversion. default is None
freq_num
number of frequencies to use in inversion
geoelectric_strike
geoelectric strike angle assuming N = 0, E = 90
masked_data
similar to data, but any masked points are now 0
mode_dict
dictionary of model modes to chose from
model_mode
model mode to use for inversion, see above
num_edi
number of stations to invert for
occam_dict
dictionary of occam parameters to use internally
occam_format
occam format of data file. default is OCCAM2MTDATA_1.0
phase_te_err
percent error in phase for TE mode. default is 5
phase_tm_err
percent error in phase for TM mode. default is 5
profile_angle
angle of profile line realtive to N = 0, E = 90
profile_line
m, b coefficients for mx+b definition of profile line
res_te_err
percent error in resistivity for TE mode. default is 10
res_tm_err
percent error in resistivity for TM mode. default is 10
error_type
‘floor’ or ‘absolute’ - option to set error as floor (i.e. maximum of the data error and a specified value) or a straight value
save_path
directory to save files to
station_list
list of station for inversion
station_locations
station locations along profile line
tipper_err
percent error in tipper. default is 5
title
title in data file.
Methods
Description
_fill_data
fills the data array that is described above
_get_data_list
gets the lines to write to data file
_get_frequencies
gets frequency list to invert for
get_profile_origin
get profile origin in UTM coordinates
mask_points
masks points in data picked from plot_mask_points
plot_mask_points
plots data responses to interactively mask data points.
plot_resonse
plots data/model responses, returns PlotResponse data type.
read_data_file
read in existing data file and fill appropriate attributes.
write_data_file
write a data file according to Data attributes
- Example Write Data File
:: >>> import mtpy.modeling.occam2d as occam2d >>> edipath = r”/home/mt/edi_files” >>> slst = [‘mt{0:03}’.format(ss) for ss in range(1, 20)] >>> ocd = occam2d.Data(edi_path=edipath, station_list=slst) >>> # model just the tm mode and tipper >>> ocd.model_mode = 3 >>> ocd.save_path = r”/home/occam/Line1/Inv1” >>> ocd.write_data_file() >>> # mask points >>> ocd.plot_mask_points() >>> ocd.mask_points()
Methods
generate_profile
()Generate linear profile by regression of station locations.
get the origin of the profile in real world coordinates
mask_from_datafile
(mask_datafn)reads a separate data file and applies mask from this data file.
mask_points
(maskpoints_obj)mask points and rewrite the data file
plot_mask_points
([data_fn, marker, ...])An interactive plotting tool to mask points an add errorbars
plot_profile
(**kwargs)Plot the projected profile line along with original station locations to make sure the line projected is correct.
plot_response
(**kwargs)plot data and model responses as apparent resistivity, phase and tipper.
project_elevation
([elevation_model])projects elevation data into the profile
read_data_file
([data_fn])Read in an existing data file and populate appropriate attributes
write_data_file
([data_fn])Write a data file.
- get_profile_origin()[source]¶
get the origin of the profile in real world coordinates
Author: Alison Kirkby (2013)
NEED TO ADAPT THIS TO THE CURRENT SETUP.
- mask_from_datafile(mask_datafn)[source]¶
reads a separate data file and applies mask from this data file. mask_datafn needs to have exactly the same frequencies, and station names must match exactly.
- mask_points(maskpoints_obj)[source]¶
mask points and rewrite the data file
NEED TO REDO THIS TO FIT THE CURRENT SETUP
- plot_mask_points(data_fn=None, marker='h', res_err_inc=0.25, phase_err_inc=0.05)[source]¶
An interactive plotting tool to mask points an add errorbars
- plot_response(**kwargs)[source]¶
plot data and model responses as apparent resistivity, phase and tipper. See PlotResponse for key words.
- class mtpy.modeling.occam2d_rewrite.Mask(edi_path=None, **kwargs)[source]¶
Allow masking of points from data file (effectively commenting them out, so the process is reversable). Inheriting from Data class.
Methods
generate_profile
()Generate linear profile by regression of station locations.
get_profile_origin
()get the origin of the profile in real world coordinates
mask_from_datafile
(mask_datafn)reads a separate data file and applies mask from this data file.
mask_points
(maskpoints_obj)mask points and rewrite the data file
plot_mask_points
([data_fn, marker, ...])An interactive plotting tool to mask points an add errorbars
plot_profile
(**kwargs)Plot the projected profile line along with original station locations to make sure the line projected is correct.
plot_response
(**kwargs)plot data and model responses as apparent resistivity, phase and tipper.
project_elevation
([elevation_model])projects elevation data into the profile
read_data_file
([data_fn])Read in an existing data file and populate appropriate attributes
write_data_file
([data_fn])Write a data file.
- class mtpy.modeling.occam2d_rewrite.Mesh(station_locations=None, **kwargs)[source]¶
deals only with the finite element mesh. Builds a finite element mesh based on given parameters defined below. The mesh reads in the station locations, finds the center and makes the relative location of the furthest left hand station 0. The mesh increases in depth logarithmically as required by the physics of MT. Also, the model extends horizontally and vertically with padding cells in order to fullfill the assumption of the forward operator that at the edges the structure is 1D. Stations are place on the horizontal nodes as required by Wannamaker’s forward operator.
Mesh has the ability to create a mesh that incorporates topography given a elevation profile. It adds more cells to the mesh with thickness z1_layer. It then sets the values of the triangular elements according to the elevation value at that location. If the elevation covers less than 50% of the triangular cell, then the cell value is set to that of air
Note
Mesh is inhereted by Regularization, so the mesh can also be be built from there, same as the example below.
Methods
add_elevation
([elevation_profile])the elevation model needs to be in relative coordinates and be a numpy.ndarray(2, num_elevation_points) where the first column is the horizontal location and the second column is the elevation at that location.
Build the finite element mesh given the parameters defined by the attributes of Mesh.
plot_mesh
(**kwargs)Plot built mesh with station locations.
read_mesh_file
(mesh_fn)reads an occam2d 2D mesh file
write_mesh_file
([save_path, basename])Write a finite element mesh file.
- add_elevation(elevation_profile=None)[source]¶
the elevation model needs to be in relative coordinates and be a numpy.ndarray(2, num_elevation_points) where the first column is the horizontal location and the second column is the elevation at that location.
If you have a elevation model use Profile to project the elevation information onto the profile line
To build the elevation I’m going to add the elevation to the top of the model which will add cells to the mesh. there might be a better way to do this, but this is the first attempt. So I’m going to assume that the first layer of the mesh without elevation is the minimum elevation and blocks will be added to max elevation at an increment according to z1_layer
Note
the elevation model should be symmetrical ie, starting at the first station and ending on the last station, so for now any elevation outside the station area will be ignored and set to the elevation of the station at the extremities. This is not ideal but works for now.
- build_mesh()[source]¶
Build the finite element mesh given the parameters defined by the attributes of Mesh. Computes relative station locations by finding the center of the station area and setting the middle to 0. Mesh blocks are built by calculating the distance between stations and putting evenly spaced blocks between the stations being close to cell_width. This places a horizontal node at the station location. If the spacing between stations is smaller than cell_width, a horizontal node is placed between the stations to be sure the model has room to change between the station.
If elevation_profile is given, add_elevation is called to add topography into the mesh.
- Populates attributes:
mesh_values
rel_station_locations
x_grid
x_nodes
z_grid
z_nodes
- Example
:: >>> import mtpy.modeling.occam2d as occcam2d >>> edipath = r”/home/mt/edi_files” >>> slist = [‘mt{0:03}’.format(ss) for ss in range(20)] >>> ocd = occam2d.Data(edi_path=edipath, station_list=slist) >>> ocd.save_path = r”/home/occam/Line1/Inv1” >>> ocd.write_data_file() >>> ocm = occam2d.Mesh(ocd.station_locations) >>> # add in elevation >>> ocm.elevation_profile = ocd.elevation_profile >>> # change number of layers >>> ocm.n_layers = 110 >>> # change cell width in station area >>> ocm.cell_width = 200 >>> ocm.build_mesh()
- plot_mesh(**kwargs)[source]¶
Plot built mesh with station locations.
Key Words
Description
depth_scale
[ ‘km’ | ‘m’ ] scale of mesh plot. default is ‘km’
fig_dpi
dots-per-inch resolution of the figure default is 300
fig_num
number of the figure instance default is ‘Mesh’
fig_size
size of figure in inches (width, height) default is [5, 5]
fs
size of font of axis tick labels, axis labels are fs+2. default is 6
ls
[ ‘-’ | ‘.’ | ‘:’ ] line style of mesh lines default is ‘-’
marker
marker of stations default is r”$lacktriangledown$”
ms
size of marker in points. default is 5
plot_triangles
[ ‘y’ | ‘n’ ] to plot mesh triangles. default is ‘n’
- class mtpy.modeling.occam2d_rewrite.Model(iter_fn=None, model_fn=None, mesh_fn=None, **kwargs)[source]¶
Read .iter file output by Occam2d. Builds the resistivity model from mesh and regularization files found from the .iter file. The resistivity model is an array(x_nodes, z_nodes) set on a regular grid, and the values of the model response are filled in according to the regularization grid. This allows for faster plotting.
Inherets Startup because they are basically the same object.
Methods
build the model from the mesh, regularization grid and model file
read_iter_file
([iter_fn])Read an iteration file.
write_iter_file
([iter_fn])write an iteration file if you need to for some reason, same as startup file
write_startup_file
([startup_fn, save_path, ...])Write a startup file based on the parameters of startup class.
- class mtpy.modeling.occam2d_rewrite.OccamPointPicker(ax_list, line_list, err_list, res_err_inc=0.05, phase_err_inc=0.02, marker='h')[source]¶
This class helps the user interactively pick points to mask and add error bars.
Methods
__call__
(event)When the function is called the mouse events will be recorder for picking points to mask or change error bars.
inAxes
(event)gets the axes that the mouse is currently in.
inFigure
(event)gets the figure number that the mouse is in
on_close
(event)close the figure with a 'q' key event and disconnect the event ids
- class mtpy.modeling.occam2d_rewrite.PlotL2(iter_fn, **kwargs)[source]¶
Plot L2 curve of iteration vs rms and rms vs roughness.
Need to only input an .iter file, will read all similar .iter files to get the rms, iteration number and roughness of all similar .iter files.
Methods
plot
()plot L2 curve
redraw plot if parameters were changed
save_figure
(save_fn[, file_format, ...])save_plot will save the figure to save_fn.
update any parameters that where changed using the built-in draw from canvas.
- 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.plotAllResponses() >>> #change line width >>> p1.lw = 2 >>> p1.redraw_plot()
- save_figure(save_fn, file_format='pdf', orientation='portrait', fig_dpi=None, close_fig='y')[source]¶
save_plot will save the figure to save_fn.
- 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.plotAllResponses() >>> [ax.grid(True, which='major') for ax in [ps1.axrte,ps1.axtep]] >>> ps1.update_plot()
- class mtpy.modeling.occam2d_rewrite.PlotMisfitPseudoSection(data_fn, resp_fn, **kwargs)[source]¶
plot a pseudo section of the data and response if given
Methods
compute misfit of MT response found from the model and the data.
plot
()plot pseudo section of data and response if given
redraw plot if parameters were changed
save_figure
(save_fn[, file_format, ...])save_plot will save the figure to save_fn.
update any parameters that where changed using the built-in draw from canvas.
- get_misfit()[source]¶
compute misfit of MT response found from the model and the data.
Need to normalize correctly
- 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.
- 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.occam2d_rewrite.PlotModel(iter_fn=None, data_fn=None, **kwargs)[source]¶
plot the 2D model found by Occam2D. The model is displayed as a meshgrid instead of model bricks. This speeds things up considerably.
Inherets the Model class to take advantage of the attributes and methods already coded.
Methods
build_model
()build the model from the mesh, regularization grid and model file
plot
()plotModel will plot the model output by occam2d in the iteration file.
read_iter_file
([iter_fn])Read an iteration file.
redraw plot if parameters were changed
save_figure
(save_fn[, file_format, ...])save_plot will save the figure to save_fn.
update any parameters that where changed using the built-in draw from canvas.
write_iter_file
([iter_fn])write an iteration file if you need to for some reason, same as startup file
write_startup_file
([startup_fn, save_path, ...])Write a startup file based on the parameters of startup class.
- plot()[source]¶
plotModel will plot the model output by occam2d in the iteration file.
- Example
>>> import mtpy.modeling.occam2d as occam2d >>> itfn = r"/home/Occam2D/Line1/Inv1/Test_15.iter" >>> model_plot = occam2d.PlotModel(itfn) >>> model_plot.ms = 20 >>> model_plot.ylimits = (0,.350) >>> model_plot.yscale = 'm' >>> model_plot.spad = .10 >>> model_plot.ypad = .125 >>> model_plot.xpad = .025 >>> model_plot.climits = (0,2.5) >>> model_plot.aspect = 'equal' >>> model_plot.redraw_plot()
- 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.plotAllResponses() >>> #change line width >>> p1.lw = 2 >>> p1.redraw_plot()
- save_figure(save_fn, file_format='pdf', orientation='portrait', fig_dpi=None, close_fig='y')[source]¶
save_plot will save the figure to save_fn.
- 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.plotAllResponses() >>> [ax.grid(True, which='major') for ax in [ps1.axrte,ps1.axtep]] >>> ps1.update_plot()
- class mtpy.modeling.occam2d_rewrite.PlotPseudoSection(data_fn, resp_fn=None, **kwargs)[source]¶
plot a pseudo section of the data and response if given
Methods
plot
()plot pseudo section of data and response if given
redraw plot if parameters were changed
save_figure
(save_fn[, file_format, ...])save_plot will save the figure to save_fn.
update any parameters that where changed using the built-in draw from canvas.
- 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.
- 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.occam2d_rewrite.PlotResponse(data_fn, resp_fn=None, **kwargs)[source]¶
Helper class to deal with plotting the MT response and occam2d model.
Methods
plot
()plot the data and model response, if given, in individual plots.
redraw plot if parameters were changed
save_figures
(save_path[, fig_fmt, fig_dpi, ...])save all the figure that are in self.fig_list
- 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')
- class mtpy.modeling.occam2d_rewrite.Profile(edi_path=None, edi_list=[], **kwargs)[source]¶
Takes data from .edi files to create a profile line for 2D modeling. Can project the stations onto a profile that is perpendicular to strike or a given profile direction.
If _rotate_to_strike is True, the impedance tensor and tipper are rotated to align with the geoelectric strike angle.
If _rotate_to_strike is True and geoelectric_strike is not given, then it is calculated using the phase tensor. First, 2D sections are estimated from the impedance tensor then the strike is estimated from the phase tensor azimuth + skew. This angle is then used to project the stations perpendicular to the strike angle.
If you want to project onto an angle not perpendicular to strike, give profile_angle and set _rotate_to_strike to False. This will project the impedance tensor and tipper to be perpendicular with the profile_angle.
Methods
Generate linear profile by regression of station locations.
plot_profile
(**kwargs)Plot the projected profile line along with original station locations to make sure the line projected is correct.
project_elevation
([elevation_model])projects elevation data into the profile
- generate_profile()[source]¶
Generate linear profile by regression of station locations.
If profile_angle is not None, then station are projected onto that line. Else, the a geoelectric strike is calculated from the data and the stations are projected onto an angle perpendicular to the estimated strike direction. If _rotate_to_strike is True, the impedance tensor and Tipper data are rotated to align with strike. Else, data is not rotated to strike.
To project stations onto a given line, set profile_angle and _rotate_to_strike to False. This will project the stations onto profile_angle and rotate the impedance tensor and tipper to be perpendicular to the profile_angle.
- plot_profile(**kwargs)[source]¶
Plot the projected profile line along with original station locations to make sure the line projected is correct.
Key Words
Description
fig_dpi
dots-per-inch resolution of figure default is 300
fig_num
number if figure instance default is ‘Projected Profile’
fig_size
size of figure in inches (width, height) default is [5, 5]
fs
[ float ] font size in points of axes tick labels axes labels are fs+2 default is 6
lc
[ string | (r, g, b) ]color of profile line (see matplotlib.line for options) default is ‘b’ – blue
lw
float, width of profile line in points default is 1
marker
[ string ] marker for stations (see matplotlib.pyplot.plot) for options
mc
[ string | (r, g, b) ] color of projected stations. default is ‘k’ – black
ms
[ float ] size of station marker default is 5
station_id
[min, max] index values for station labels default is None
- Example
:: >>> edipath = r”/home/mt/edi_files” >>> pr = occam2d.Profile(edi_path=edipath) >>> pr.generate_profile() >>> # set station labels to only be from 1st to 4th index >>> # of station name >>> pr.plot_profile(station_id=[0,4])
- class mtpy.modeling.occam2d_rewrite.Regularization(station_locations=None, **kwargs)[source]¶
Creates a regularization grid based on Mesh. Note that Mesh is inherited by Regularization, therefore the intended use is to build a mesh with the Regularization class.
The regularization grid is what Occam calculates the inverse model on. Setup is tricky and can be painful, as you can see it is not quite fully functional yet, as it cannot incorporate topography yet. It seems like you’d like to have the regularization setup so that your target depth is covered well, in that the regularization blocks to this depth are sufficiently small to resolve resistivity structure at that depth. Finally, you want the regularization to go to a half space at the bottom, basically one giant block.
Methods
add_elevation
([elevation_profile])the elevation model needs to be in relative coordinates and be a numpy.ndarray(2, num_elevation_points) where the first column is the horizontal location and the second column is the elevation at that location.
build_mesh
()Build the finite element mesh given the parameters defined by the attributes of Mesh.
Builds larger boxes around existing mesh blocks for the regularization.
estimate the number of free parameters in model mesh.
plot_mesh
(**kwargs)Plot built mesh with station locations.
read_mesh_file
(mesh_fn)reads an occam2d 2D mesh file
read_regularization_file
(reg_fn)Read in a regularization file and populate attributes:
write_mesh_file
([save_path, basename])Write a finite element mesh file.
write_regularization_file
([reg_fn, ...])Write a regularization file for input into occam.
- build_regularization()[source]¶
Builds larger boxes around existing mesh blocks for the regularization. As the model deepens the regularization boxes get larger.
The regularization boxes are merged mesh cells as prescribed by the Occam method.
- get_num_free_params()[source]¶
estimate the number of free parameters in model mesh.
I’m assuming that if there are any fixed parameters in the block, then that model block is assumed to be fixed. Not sure if this is right cause there is no documentation.
DOES NOT WORK YET
- read_regularization_file(reg_fn)[source]¶
- Read in a regularization file and populate attributes:
binding_offset
mesh_fn
model_columns
model_rows
prejudice_fn
statics_fn
- write_regularization_file(reg_fn=None, reg_basename=None, statics_fn='none', prejudice_fn='none', save_path=None)[source]¶
Write a regularization file for input into occam.
Calls build_regularization if build_regularization has not already been called.
if reg_fn is None, then file is written to save_path/reg_basename
- class mtpy.modeling.occam2d_rewrite.Response(resp_fn=None, **kwargs)[source]¶
Reads .resp file output by Occam. Similar structure to Data.data.
If resp_fn is given in the initialization of Response, read_response_file is called.
Methods
read_response_file
([resp_fn])read in response file and put into a list of dictionaries similar to Data
- class mtpy.modeling.occam2d_rewrite.Run[source]¶
Run Occam2D by system call.
Future plan: implement Occam in Python and call it from here directly.
- class mtpy.modeling.occam2d_rewrite.Startup(**kwargs)[source]¶
Reads and writes the startup file for Occam2D.
Note
Be sure to look at the Occam 2D documentation for description of all parameters
Key Words/Attributes
Description
data_fn
full path to data file
date_time
date and time the startup file was written
debug_level
[ 0 | 1 | 2 ] see occam documentation default is 1
description
brief description of inversion run default is ‘startup created by mtpy’
diagonal_penalties
penalties on diagonal terms default is 0
format
Occam file format default is ‘OCCAMITER_FLEX’
iteration
current iteration number default is 0
iterations_to_run
maximum number of iterations to run default is 20
lagrange_value
starting lagrange value default is 5
misfit_reached
[ 0 | 1 ] 0 if misfit has been reached, 1 if it has. default is 0
misfit_value
current misfit value. default is 1000
model_fn
full path to model file
model_limits
limits on model resistivity values default is None
model_value_steps
limits on the step size of model values default is None
model_values
np.ndarray(num_free_params) of model values
param_count
number of free parameters in model
resistivity_start
starting resistivity value. If model_values is not given, then all values with in model_values array will be set to resistivity_start
roughness_type
[ 0 | 1 | 2 ] type of roughness default is 1
roughness_value
current roughness value. default is 1E10
save_path
directory path to save startup file to default is current working directory
startup_basename
basename of startup file name. default is Occam2DStartup
startup_fn
full path to startup file. default is save_path/startup_basename
stepsize_count
max number of iterations per step default is 8
target_misfit
target misfit value. default is 1.
- Example
>>> startup = occam2d.Startup() >>> startup.data_fn = ocd.data_fn >>> startup.model_fn = profile.reg_fn >>> startup.param_count = profile.num_free_params >>> startup.save_path = r"/home/occam2d/Line1/Inv1"
Methods
write_startup_file
([startup_fn, save_path, ...])Write a startup file based on the parameters of startup class.
Module Winglink¶
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]¶
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.
Methods
compute misfit of MT response found from the model and the data.
plot
()plot pseudo section of data and response if given
redraw plot if parameters were changed
save_figure
(save_fn[, file_format, ...])save_plot will save the figure to save_fn.
update any parameters that where changed using the built-in draw from canvas.
- get_misfit()[source]¶
compute misfit of MT response found from the model and the data.
Need to normalize correctly
- 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.
- 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]¶
plot a pseudo section of the data and response if given
Methods
plot
()plot pseudo section of data and response if given
redraw plot if parameters were changed
save_figure
(save_fn[, file_format, ...])save_plot will save the figure to save_fn.
update any parameters that where changed using the built-in draw from canvas.
- 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.
- 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]¶
Helper class to deal with plotting the MT response and occam2d model.
Methods
plot
()plot the data and model response, if given, in individual plots.
redraw plot if parameters were changed
save_figures
(save_path[, fig_fmt, fig_dpi, ...])save all the figure that are in self.fig_list
- 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')
Module WS3DINV¶
- Deals with input and output files for ws3dinv written by:
Siripunvaraporn, W.; Egbert, G.; Lenbury, Y. & Uyeshima, M. Three-dimensional magnetotelluric inversion: data-space method Physics of The Earth and Planetary Interiors, 2005, 150, 3-14 * Dependencies: matplotlib 1.3.x, numpy 1.7.x, scipy 0.13
and evtk if vtk files want to be written.
The intended use or workflow is something like this for getting started:
- Making input files
>>> import mtpy.modeling.ws3dinv as ws >>> 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 edi_path >>> ... if edi.find('.edi') > 0] >>> #2) make a grid from the stations themselves with 200m cell spacing >>> wsmesh = ws.WSMesh(edi_list=edi_list, cell_size_east=200, >>> ... cell_size_north=200) >>> wsmesh.make_mesh() >>> # check to see if the mesh is what you think it should be >>> wsmesh.plot_mesh() >>> # all is good write the mesh file >>> wsmesh.write_initial_file(save_path=r"/home/ws3dinv/Inv1") >>> # note this will write a file with relative station locations >>> #change the starting model to be different than a halfspace >>> mm = ws.WS3DModelManipulator(initial_fn=wsmesh.initial_fn) >>> # an interactive gui will pop up to change the resistivity model >>> #once finished write a new initial file >>> mm.rewrite_initial_file() >>> #3) write data file >>> wsdata = ws.WSData(edi_list=edi_list, station_fn=wsmesh.station_fn) >>> wsdata.write_data_file() >>> #4) plot mt response to make sure everything looks ok >>> rp = ws.PlotResponse(data_fn=wsdata.data_fn) >>> #5) make startup file >>> sws = ws.WSStartup(data_fn=wsdata.data_fn, initial_fn=mm.new_initial_fn)
- checking the model and response
>>> mfn = r"/home/ws3dinv/Inv1/test_model.01" >>> dfn = r"/home/ws3dinv/Inv1/WSDataFile.dat" >>> rfn = r"/home/ws3dinv/Inv1/test_resp.01" >>> sfn = r"/home/ws3dinv/Inv1/WS_Sation_Locations.txt" >>> # plot the data vs. model response >>> rp = ws.PlotResponse(data_fn=dfn, resp_fn=rfn, station_fn=sfn) >>> # plot model slices where you can interactively step through >>> ds = ws.PlotSlices(model_fn=mfn, station_fn=sfn) >>> # plot phase tensor ellipses on top of depth slices >>> ptm = ws.PlotPTMaps(data_fn=dfn, resp_fn=rfn, model_fn=mfn) >>> #write files for 3D visualization in Paraview or Mayavi >>> ws.write_vtk_files(mfn, sfn, r"/home/ParaviewFiles")
Created on Sun Aug 25 18:41:15 2013
@author: jpeacock-pr
- class mtpy.modeling.ws3dinv.PlotDepthSlice(model_fn=None, data_fn=None, station_fn=None, initial_fn=None, **kwargs)[source]¶
Plots depth slices of resistivity model
- Example
>>> import mtpy.modeling.ws3dinv as ws >>> mfn = r"/home/MT/ws3dinv/Inv1/Test_model.00" >>> sfn = r"/home/MT/ws3dinv/Inv1/WSStationLocations.txt" >>> # plot just first layer to check the formating >>> pds = ws.PlotDepthSlice(model_fn=mfn, station_fn=sfn, >>> ... depth_index=0, save_plots='n') >>> #move color bar up >>> pds.cb_location >>> (0.64500000000000002, 0.14999999999999997, 0.3, 0.025) >>> pds.cb_location = (.645, .175, .3, .025) >>> pds.redraw_plot() >>> #looks good now plot all depth slices and save them to a folder >>> pds.save_path = r"/home/MT/ws3dinv/Inv1/DepthSlices" >>> pds.depth_index = None >>> pds.save_plots = 'y' >>> pds.redraw_plot()
Attributes
Description
cb_location
location of color bar (x, y, width, height) default is None, automatically locates
cb_orientation
[ ‘vertical’ | ‘horizontal’ ] default is horizontal
cb_pad
padding between axes and colorbar default is None
cb_shrink
percentage to shrink colorbar by default is None
climits
(min, max) of resistivity color on log scale default is (0, 4)
cmap
name of color map default is ‘jet_r’
data_fn
full path to data file
depth_index
integer value of depth slice index, shallowest layer is 0
dscale
scaling parameter depending on map_scale
ew_limits
(min, max) plot limits in e-w direction in map_scale units. default is None, sets viewing area to the station area
fig_aspect
aspect ratio of plot. default is 1
fig_dpi
resolution of figure in dots-per-inch. default is 300
fig_list
list of matplotlib.figure instances for each depth slice
fig_size
[width, height] in inches of figure size default is [6, 6]
font_size
size of ticklabel font in points, labels are font_size+2. default is 7
grid_east
relative location of grid nodes in e-w direction in map_scale units
grid_north
relative location of grid nodes in n-s direction in map_scale units
grid_z
relative location of grid nodes in z direction in map_scale units
initial_fn
full path to initial file
map_scale
[ ‘km’ | ‘m’ ] distance units of map. default is km
mesh_east
np.meshgrid(grid_east, grid_north, indexing=’ij’)
mesh_north
np.meshgrid(grid_east, grid_north, indexing=’ij’)
model_fn
full path to model file
nodes_east
relative distance betwen nodes in e-w direction in map_scale units
nodes_north
relative distance betwen nodes in n-s direction in map_scale units
nodes_z
relative distance betwen nodes in z direction in map_scale units
ns_limits
(min, max) plot limits in n-s direction in map_scale units. default is None, sets viewing area to the station area
plot_grid
[ ‘y’ | ‘n’ ] ‘y’ to plot mesh grid lines. default is ‘n’
plot_yn
[ ‘y’ | ‘n’ ] ‘y’ to plot on instantiation
res_model
np.ndarray(n_north, n_east, n_vertical) of model resistivity values in linear scale
save_path
path to save figures to
save_plots
[ ‘y’ | ‘n’ ] ‘y’ to save depth slices to save_path
station_east
location of stations in east direction in map_scale units
station_fn
full path to station locations file
station_names
station names
station_north
location of station in north direction in map_scale units
subplot_bottom
distance between axes and bottom of figure window
subplot_left
distance between axes and left of figure window
subplot_right
distance between axes and right of figure window
subplot_top
distance between axes and top of figure window
title
titiel of plot default is depth of slice
xminorticks
location of xminorticks
yminorticks
location of yminorticks
Methods
plot
()plot depth slices
read in the files to get appropriate information
redraw plot if parameters were changed
update_plot
(fig)update any parameters that where changed using the built-in draw from canvas.
- 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.plotAllResponses() >>> #change line width >>> p1.lw = 2 >>> p1.redraw_plot()
- update_plot(fig)[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.plotAllResponses() >>> [ax.grid(True, which='major') for ax in [ps1.axrte,ps1.axtep]] >>> ps1.update_plot()
- class mtpy.modeling.ws3dinv.PlotPTMaps(data_fn=None, resp_fn=None, station_fn=None, model_fn=None, initial_fn=None, **kwargs)[source]¶
Plot phase tensor maps including residual pt if response file is input.
- Plot only data for one period
>>> import mtpy.modeling.ws3dinv as ws >>> dfn = r"/home/MT/ws3dinv/Inv1/WSDataFile.dat" >>> ptm = ws.PlotPTMaps(data_fn=dfn, plot_period_list=[0])
- Plot data and model response
>>> import mtpy.modeling.ws3dinv as ws >>> dfn = r"/home/MT/ws3dinv/Inv1/WSDataFile.dat" >>> rfn = r"/home/MT/ws3dinv/Inv1/Test_resp.00" >>> mfn = r"/home/MT/ws3dinv/Inv1/Test_model.00" >>> ptm = ws.PlotPTMaps(data_fn=dfn, resp_fn=rfn, model_fn=mfn, >>> ... plot_period_list=[0]) >>> # adjust colorbar >>> ptm.cb_res_pad = 1.25 >>> ptm.redraw_plot()
Attributes
Description
cb_pt_pad
percentage from top of axes to place pt color bar. default is .90
cb_res_pad
percentage from bottom of axes to place resistivity color bar. default is 1.2
cb_residual_tick_step
tick step for residual pt. default is 3
cb_tick_step
tick step for phase tensor color bar, default is 45
data
np.ndarray(n_station, n_periods, 2, 2) impedance tensors for station data
data_fn
full path to data fle
dscale
scaling parameter depending on map_scale
ellipse_cmap
color map for pt ellipses. default is mt_bl2gr2rd
ellipse_colorby
- [ ‘skew’ | ‘skew_seg’ | ‘phimin’ | ‘phimax’|
‘phidet’ | ‘ellipticity’ ] parameter to color ellipses by. default is ‘phimin’
ellipse_range
(min, max, step) min and max of colormap, need to input step if plotting skew_seg
ellipse_size
relative size of ellipses in map_scale
ew_limits
limits of plot in e-w direction in map_scale units. default is None, scales to station area
fig_aspect
aspect of figure. default is 1
fig_dpi
resolution in dots-per-inch. default is 300
fig_list
list of matplotlib.figure instances for each figure plotted.
fig_size
[width, height] in inches of figure window default is [6, 6]
font_size
font size of ticklabels, axes labels are font_size+2. default is 7
grid_east
relative location of grid nodes in e-w direction in map_scale units
grid_north
relative location of grid nodes in n-s direction in map_scale units
grid_z
relative location of grid nodes in z direction in map_scale units
initial_fn
full path to initial file
map_scale
[ ‘km’ | ‘m’ ] distance units of map. default is km
mesh_east
np.meshgrid(grid_east, grid_north, indexing=’ij’)
mesh_north
np.meshgrid(grid_east, grid_north, indexing=’ij’)
model_fn
full path to model file
nodes_east
relative distance betwen nodes in e-w direction in map_scale units
nodes_north
relative distance betwen nodes in n-s direction in map_scale units
nodes_z
relative distance betwen nodes in z direction in map_scale units
ns_limits
(min, max) limits of plot in n-s direction default is None, viewing area is station area
pad_east
padding from extreme stations in east direction
pad_north
padding from extreme stations in north direction
period_list
list of periods from data
plot_grid
[ ‘y’ | ‘n’ ] ‘y’ to plot grid lines default is ‘n’
plot_period_list
list of period index values to plot default is None
plot_yn
[‘y’ | ‘n’ ] ‘y’ to plot on instantiation default is ‘y’
res_cmap
colormap for resisitivity values. default is ‘jet_r’
res_limits
(min, max) resistivity limits in log scale default is (0, 4)
res_model
np.ndarray(n_north, n_east, n_vertical) of model resistivity values in linear scale
residual_cmap
color map for pt residuals. default is ‘mt_wh2or’
resp
np.ndarray(n_stations, n_periods, 2, 2) impedance tensors for model response
resp_fn
full path to response file
save_path
directory to save figures to
save_plots
[ ‘y’ | ‘n’ ] ‘y’ to save plots to save_path
station_east
location of stations in east direction in map_scale units
station_fn
full path to station locations file
station_names
station names
station_north
location of station in north direction in map_scale units
subplot_bottom
distance between axes and bottom of figure window
subplot_left
distance between axes and left of figure window
subplot_right
distance between axes and right of figure window
subplot_top
distance between axes and top of figure window
title
titiel of plot default is depth of slice
xminorticks
location of xminorticks
yminorticks
location of yminorticks
Methods
plot
()plot phase tensor maps for data and or response, each figure is of a different period.
redraw plot if parameters were changed
save_figure
([save_path, fig_dpi, ...])save_figure will save the figure to save_fn.
- plot()[source]¶
plot phase tensor maps for data and or response, each figure is of a different period. If response is input a third column is added which is the residual phase tensor showing where the model is not fitting the data well. The data is plotted in km in units of ohm-m.
- Inputs:
data_fn = full path to data file resp_fn = full path to response file, if none just plots data sites_fn = full path to sites file periodlst = indicies of periods you want to plot esize = size of ellipses as:
0 = phase tensor ellipse 1 = phase tensor residual 2 = resistivity tensor ellipse 3 = resistivity tensor residual
ecolor = ‘phimin’ for coloring with phimin or ‘beta’ for beta coloring colormm = list of min and max coloring for plot, list as follows:
0 = phase tensor min and max for ecolor in degrees 1 = phase tensor residual min and max [0,1] 2 = resistivity tensor coloring as resistivity on log scale 3 = resistivity tensor residual coloring as resistivity on
linear scale
xpad = padding of map from stations at extremities (km) units = ‘mv’ to convert to Ohm-m dpi = dots per inch of figure
- 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.plotAllResponses() >>> #change line width >>> p1.lw = 2 >>> p1.redraw_plot()
- class mtpy.modeling.ws3dinv.PlotResponse(data_fn=None, resp_fn=None, station_fn=None, **kwargs)[source]¶
plot data and response
- Example
>>> import mtpy.modeling.ws3dinv as ws >>> dfn = r"/home/MT/ws3dinv/Inv1/WSDataFile.dat" >>> rfn = r"/home/MT/ws3dinv/Inv1/Test_resp.00" >>> sfn = r"/home/MT/ws3dinv/Inv1/WSStationLocations.txt" >>> wsrp = ws.PlotResponse(data_fn=dfn, resp_fn=rfn, station_fn=sfn) >>> # plot only the TE and TM modes >>> wsrp.plot_component = 2 >>> wsrp.redraw_plot()
Attributes
Description
color_mode
[ ‘color’ | ‘bw’ ] color or black and white plots
cted
color for data TE mode
ctem
color for data TM mode
ctmd
color for model TE mode
ctmm
color for model TM mode
data_fn
full path to data file
data_object
WSResponse instance
e_capsize
cap size of error bars in points (default is .5)
e_capthick
cap thickness of error bars in points (default is 1)
fig_dpi
resolution of figure in dots-per-inch (300)
fig_list
list of matplotlib.figure instances for plots
fig_size
size of figure in inches (default is [6, 6])
font_size
size of font for tick labels, axes labels are font_size+2 (default is 7)
legend_border_axes_pad
padding between legend box and axes
legend_border_pad
padding between border of legend and symbols
legend_handle_text_pad
padding between text labels and symbols of legend
legend_label_spacing
padding between labels
legend_loc
location of legend
legend_marker_scale
scale of symbols in legend
lw
line width response curves (default is .5)
ms
size of markers (default is 1.5)
mted
marker for data TE mode
mtem
marker for data TM mode
mtmd
marker for model TE mode
mtmm
marker for model TM mode
phase_limits
limits of phase
plot_component
[ 2 | 4 ] 2 for TE and TM or 4 for all components
plot_style
[ 1 | 2 ] 1 to plot each mode in a seperate subplot and 2 to plot xx, xy and yx, yy in same plots
plot_type
[ ‘1’ | list of station name ] ‘1’ to plot all stations in data file or input a list of station names to plot if station_fn is input, otherwise input a list of integers associated with the index with in the data file, ie 2 for 2nd station
plot_z
[ True | False ] default is True to plot impedance, False for plotting resistivity and phase
plot_yn
[ ‘n’ | ‘y’ ] to plot on instantiation
res_limits
limits of resistivity in linear scale
resp_fn
full path to response file
resp_object
WSResponse object for resp_fn, or list of WSResponse objects if resp_fn is a list of response files
station_fn
full path to station file written by WSStation
subplot_bottom
space between axes and bottom of figure
subplot_hspace
space between subplots in vertical direction
subplot_left
space between axes and left of figure
subplot_right
space between axes and right of figure
subplot_top
space between axes and top of figure
subplot_wspace
space between subplots in horizontal direction
Methods
plot
()plot_errorbar
(ax, period, data, error, ...)convinience function to make an error bar instance
redraw plot if parameters were changed
save_figure
(save_fn[, file_format, ...])save_plot will save the figure to save_fn.
update any parameters that where changed using the built-in draw from canvas.
- plot_errorbar(ax, period, data, error, color, marker)[source]¶
convinience function to make an error bar instance
- 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.plotAllResponses() >>> #change line width >>> p1.lw = 2 >>> p1.redraw_plot()
- save_figure(save_fn, file_format='pdf', orientation='portrait', fig_dpi=None, close_fig='y')[source]¶
save_plot will save the figure to save_fn.
- 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.plotAllResponses() >>> [ax.grid(True, which='major') for ax in [ps1.axrte,ps1.axtep]] >>> ps1.update_plot()
- class mtpy.modeling.ws3dinv.PlotSlices(model_fn, data_fn=None, station_fn=None, initial_fn=None, **kwargs)[source]¶
plot all slices and be able to scroll through the model
- Example
>>> import mtpy.modeling.ws3dinv as ws >>> mfn = r"/home/MT/ws3dinv/Inv1/Test_model.00" >>> sfn = r"/home/MT/ws3dinv/Inv1/WSStationLocations.txt" >>> # plot just first layer to check the formating >>> pds = ws.PlotSlices(model_fn=mfn, station_fn=sfn)
Buttons
Description
‘e’
moves n-s slice east by one model block
‘w’
moves n-s slice west by one model block
‘n’
moves e-w slice north by one model block
‘m’
moves e-w slice south by one model block
‘d’
moves depth slice down by one model block
‘u’
moves depth slice up by one model block
Attributes
Description
ax_en
matplotlib.axes instance for depth slice map view
ax_ez
matplotlib.axes instance for e-w slice
ax_map
matplotlib.axes instance for location map
ax_nz
matplotlib.axes instance for n-s slice
climits
(min , max) color limits on resistivity in log scale. default is (0, 4)
cmap
name of color map for resisitiviy. default is ‘jet_r’
data_fn
full path to data file name
dscale
scaling parameter depending on map_scale
east_line_xlist
list of line nodes of east grid for faster plotting
east_line_ylist
list of line nodes of east grid for faster plotting
ew_limits
(min, max) limits of e-w in map_scale units default is None and scales to station area
fig
matplotlib.figure instance for figure
fig_aspect
aspect ratio of plots. default is 1
fig_dpi
resolution of figure in dots-per-inch default is 300
fig_num
figure instance number
fig_size
[width, height] of figure window. default is [6,6]
font_dict
dictionary of font keywords, internally created
font_size
size of ticklables in points, axes labes are font_size+2. default is 7
grid_east
relative location of grid nodes in e-w direction in map_scale units
grid_north
relative location of grid nodes in n-s direction in map_scale units
grid_z
relative location of grid nodes in z direction in map_scale units
index_east
index value of grid_east being plotted
index_north
index value of grid_north being plotted
index_vertical
index value of grid_z being plotted
initial_fn
full path to initial file
key_press
matplotlib.canvas.connect instance
map_scale
[ ‘m’ | ‘km’ ] scale of map. default is km
mesh_east
np.meshgrid(grid_east, grid_north)[0]
mesh_en_east
np.meshgrid(grid_east, grid_north)[0]
mesh_en_north
np.meshgrid(grid_east, grid_north)[1]
mesh_ez_east
np.meshgrid(grid_east, grid_z)[0]
mesh_ez_vertical
np.meshgrid(grid_east, grid_z)[1]
mesh_north
np.meshgrid(grid_east, grid_north)[1]
mesh_nz_north
np.meshgrid(grid_north, grid_z)[0]
mesh_nz_vertical
np.meshgrid(grid_north, grid_z)[1]
model_fn
full path to model file
ms
size of station markers in points. default is 2
nodes_east
relative distance betwen nodes in e-w direction in map_scale units
nodes_north
relative distance betwen nodes in n-s direction in map_scale units
nodes_z
relative distance betwen nodes in z direction in map_scale units
north_line_xlist
list of line nodes north grid for faster plotting
north_line_ylist
list of line nodes north grid for faster plotting
ns_limits
(min, max) limits of plots in n-s direction default is None, set veiwing area to station area
plot_yn
[ ‘y’ | ‘n’ ] ‘y’ to plot on instantiation default is ‘y’
res_model
np.ndarray(n_north, n_east, n_vertical) of model resistivity values in linear scale
station_color
color of station marker. default is black
station_dict_east
location of stations for each east grid row
station_dict_north
location of stations for each north grid row
station_east
location of stations in east direction
station_fn
full path to station file
station_font_color
color of station label
station_font_pad
padding between station marker and label
station_font_rotation
angle of station label
station_font_size
font size of station label
station_font_weight
weight of font for station label
station_id
[min, max] index values for station labels
station_marker
station marker
station_names
name of stations
station_north
location of stations in north direction
subplot_bottom
distance between axes and bottom of figure window
subplot_hspace
distance between subplots in vertical direction
subplot_left
distance between axes and left of figure window
subplot_right
distance between axes and right of figure window
subplot_top
distance between axes and top of figure window
subplot_wspace
distance between subplots in horizontal direction
title
title of plot
z_limits
(min, max) limits in vertical direction,
Methods
get the grid line on which a station resides for plotting
on_key_press
(event)on a key press change the slices
plot
()plot:
read in the files to get appropriate information
redraw plot if parameters were changed
save_figure
([save_fn, fig_dpi, file_format, ...])save_figure will save the figure to save_fn.
- 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.plotAllResponses() >>> #change line width >>> p1.lw = 2 >>> p1.redraw_plot()
- class mtpy.modeling.ws3dinv.WSData(**kwargs)[source]¶
Includes tools for reading and writing data files intended to be used with ws3dinv.
- Example
>>> import mtpy.modeling.ws3dinv as ws >>> import os >>> edi_path = r"/home/EDI_Files" >>> edi_list = [os.path.join(edi_path, edi) for edi in edi_path >>> ... if edi.find('.edi') > 0] >>> # create an evenly space period list in log space >>> p_list = np.logspace(np.log10(.001), np.log10(1000), 12) >>> wsdata = ws.WSData(edi_list=edi_list, period_list=p_list, >>> ... station_fn=r"/home/stations.txt") >>> wsdata.write_data_file()
Attributes
Description
data
- numpy structured array with keys:
data_fn
full path to data file
edi_list
list of edi files used to make data file
n_z
[ 4 | 8 ] number of impedance tensor elements default is 8
ncol
number of columns in out file from winglink default is 5
period_list
list of periods to invert for
ptol
if periods in edi files don’t match period_list then program looks for periods within ptol defualt is .15 or 15 percent
rotation_angle
Angle to rotate the data relative to north. Here the angle is measure clockwise from North, Assuming North is 0 and East is 90. Rotating data, and grid to align with regional geoelectric strike can improve the inversion. default is None
save_path
path to save the data file
station_fn
full path to station file written by WSStation
station_locations
- numpy structured array for station locations keys:
station –> station name
- east –> relative eastern location in
grid
- north –> relative northern location in
grid
if input a station file is written
station_east
relative locations of station in east direction
station_north
relative locations of station in north direction
station_names
names of stations
units
[ ‘mv’ | ‘else’ ] units of Z, needs to be mv for ws3dinv. default is ‘mv’
wl_out_fn
Winglink .out file which describes a 3D grid
wl_site_fn
Wingling .sites file which gives station locations
z_data
impedance tensors of data with shape: (n_station, n_periods, 2, 2)
z_data_err
error of data impedance tensors with error map applied, shape (n_stations, n_periods, 2, 2)
z_err
[ float | ‘data’ ] ‘data’ to set errors as data errors or give a percent error to impedance tensor elements default is .05 or 5% if given as percent, ie. 5% then it is converted to .05.
z_err_floor
percent error floor, anything below this error will be set to z_err_floor. default is None
z_err_map
[zxx, zxy, zyx, zyy] for n_z = 8 [zxy, zyx] for n_z = 4 Value in percent to multiply the error by, which give the user power to down weight bad data, so the resulting error will be z_err_map*z_err
Methods
Description
build_data
builds the data from .edi files
write_data_file
writes a data file from attribute data. This way you can read in a data file, change some parameters and rewrite.
read_data_file
reads in a ws3dinv data file
Methods
Builds the data from .edi files to be written into a data file
compute the errors from the given attributes
read_data_file
([data_fn, wl_sites_fn, ...])read in data file
write_data_file
(**kwargs)Writes a data file based on the attribute data
- class mtpy.modeling.ws3dinv.WSMesh(edi_list=None, **kwargs)[source]¶
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.
- Example
>>> import mtpy.modeling.ws3dinv as ws >>> 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 edi_path >>> ... if edi.find('.edi') > 0] >>> #2) make a grid from the stations themselves with 200m cell spacing >>> wsmesh = ws.WSMesh(edi_list=edi_list, cell_size_east=200, >>> ... cell_size_north=200) >>> wsmesh.make_mesh() >>> # check to see if the mesh is what you think it should be >>> wsmesh.plot_mesh() >>> # all is good write the mesh file >>> wsmesh.write_initial_file(save_path=r"/home/ws3dinv/Inv1")
Attributes
Description
cell_size_east
mesh block width in east direction default is 500
cell_size_north
mesh block width in north direction default is 500
edi_list
list of .edi files to invert for
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
initial_fn
full path to initial file name
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 5
pad_north
number of cells for padding on S and N sides default is 5
pad_root_east
padding cells E & W will be pad_root_east**(x)
pad_root_north
padding cells N & S will be pad_root_north**(x)
pad_z
number of cells for padding at bottom default is 5
res_list
list of resistivity values for starting model
res_model
starting resistivity model
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
station_fn
full path to station file
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
Methods
Description
make_mesh
makes a mesh from the given specifications
plot_mesh
plots mesh to make sure everything is good
write_initial_file
writes an initial model file that includes the mesh
Methods
convert the resistivity model that is in ohm-m to integer values corresponding to res_list
create finite element mesh according to parameters set.
plot_mesh
([east_limits, north_limits, z_limits])read_initial_file
(initial_fn)read an initial file and return the pertinent information including grid positions in coordinates relative to the center point (0,0) and starting model.
write_initial_file
(**kwargs)will write an initial file for wsinv3d.
- convert_model_to_int()[source]¶
convert the resistivity model that is in ohm-m to integer values corresponding to res_list
- make_mesh()[source]¶
create finite element mesh according to parameters set.
The mesh is built by first finding the center of the station area. Then cells are added in the north and east direction with width cell_size_east and cell_size_north to the extremeties of the station area. Padding cells are then added to extend the model to reduce edge effects. The number of cells are pad_east and pad_north and the increase in size is by pad_root_east and pad_root_north. The station locations are then computed as the center of the nearest cell as required by the code.
The vertical cells are built to increase in size exponentially with depth. The first cell depth is first_layer_thickness and should be about 1/10th the shortest skin depth. The layers then increase on a log scale to z_target_depth. Then the model is padded with pad_z number of cells to extend the depth of the model.
- padding = np.round(cell_size_east*pad_root_east**np.arange(start=.5,
stop=3, step=3./pad_east))+west
- read_initial_file(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.
- write_initial_file(**kwargs)[source]¶
will write an initial file for wsinv3d.
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.
- class mtpy.modeling.ws3dinv.WSModel(model_fn=None)[source]¶
Reads in model file and fills necessary attributes.
- Example
>>> mfn = r"/home/ws3dinv/test_model.00" >>> wsmodel = ws.WSModel(mfn) >>> wsmodel.write_vtk_file(r"/home/ParaviewFiles")
Attributes
Description
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
iteration_number
iteration number of the inversion
lagrange
lagrange multiplier
model_fn
full path to model file
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
res_model
starting resistivity model
rms
root mean squared error of data and model
Methods
Description
read_model_file
read model file and fill attributes
write_vtk_file
write a vtk structured grid file for resistivity model
Methods
read in a model file as x-north, y-east, z-positive down
write_vtk_file
- class mtpy.modeling.ws3dinv.WSModelManipulator(model_fn=None, initial_fn=None, data_fn=None, **kwargs)[source]¶
will plot a model from wsinv3d or init file so the user can manipulate the resistivity values relatively easily. At the moment only plotted in map view.
- Example
:: >>> import mtpy.modeling.ws3dinv as ws >>> initial_fn = r”/home/MT/ws3dinv/Inv1/WSInitialFile” >>> mm = ws.WSModelManipulator(initial_fn=initial_fn)
Buttons
Description
‘=’
increase depth to next vertical node (deeper)
‘-’
decrease depth to next vertical node (shallower)
‘q’
quit the plot, rewrites initial file when pressed
‘a’
copies the above horizontal layer to the present layer
‘b’
copies the below horizonal layer to present layer
‘u’
undo previous change
Attributes
Description
ax1
matplotlib.axes instance for mesh plot of the model
ax2
matplotlib.axes instance of colorbar
cb
matplotlib.colorbar instance for colorbar
cid_depth
matplotlib.canvas.connect for depth
cmap
matplotlib.colormap instance
cmax
maximum value of resistivity for colorbar. (linear)
cmin
minimum value of resistivity for colorbar (linear)
data_fn
full path fo data file
depth_index
integer value of depth slice for plotting
dpi
resolution of figure in dots-per-inch
dscale
depth scaling, computed internally
east_line_xlist
list of east mesh lines for faster plotting
east_line_ylist
list of east mesh lines for faster plotting
fdict
dictionary of font properties
fig
matplotlib.figure instance
fig_num
number of figure instance
fig_size
size of figure in inches
font_size
size of font in points
grid_east
location of east nodes in relative coordinates
grid_north
location of north nodes in relative coordinates
grid_z
location of vertical nodes in relative coordinates
initial_fn
full path to initial file
m_height
mean height of horizontal cells
m_width
mean width of horizontal cells
map_scale
[ ‘m’ | ‘km’ ] scale of map
mesh_east
np.meshgrid of east, north
mesh_north
np.meshgrid of east, north
mesh_plot
matplotlib.axes.pcolormesh instance
model_fn
full path to model file
new_initial_fn
full path to new initial file
nodes_east
spacing between east nodes
nodes_north
spacing between north nodes
nodes_z
spacing between vertical nodes
north_line_xlist
list of coordinates of north nodes for faster plotting
north_line_ylist
list of coordinates of north nodes for faster plotting
plot_yn
[ ‘y’ | ‘n’ ] plot on instantiation
radio_res
matplotlib.widget.radio instance for change resistivity
rect_selector
matplotlib.widget.rect_selector
res
np.ndarray(nx, ny, nz) for model in linear resistivity
res_copy
copy of res for undo
res_dict
dictionary of segmented resistivity values
res_list
list of resistivity values for model linear scale
res_model
np.ndarray(nx, ny, nz) of resistivity values from res_list (linear scale)
res_model_int
np.ndarray(nx, ny, nz) of integer values corresponding to res_list for initial model
res_value
current resistivty value of radio_res
save_path
path to save initial file to
station_east
station locations in east direction
station_north
station locations in north direction
xlimits
limits of plot in e-w direction
ylimits
limits of plot in n-s direction
Methods
change_model_res
(xchange, ychange)change resistivity values of resistivity model
convert the resistivity model that is in ohm-m to integer values corresponding to res_list
convert_res_to_model
(res_array)converts an output model into an array of segmented valued according to res_list.
plot
()plots the model with:
reads in initial file or model file and set attributes:
rect_onselect
(eclick, erelease)on selecting a rectangle change the colors to the resistivity values
redraws the plot
rewrite_initial_file
([save_path])write an initial file for wsinv3d from the model created.
set_res_list
(res_list)on setting res_list also set the res_dict to correspond
set_res_value
- convert_model_to_int()[source]¶
convert the resistivity model that is in ohm-m to integer values corresponding to res_list
- convert_res_to_model(res_array)[source]¶
converts an output model into an array of segmented valued according to res_list.
output is an array of segemented resistivity values in ohm-m (linear)
- plot()[source]¶
- plots the model with:
-a radio dial for depth slice -radio dial for resistivity value
- read_file()[source]¶
- reads in initial file or model file and set attributes:
-resmodel -northrid -eastrid -zgrid -res_list if initial file
- rect_onselect(eclick, erelease)[source]¶
on selecting a rectangle change the colors to the resistivity values
- class mtpy.modeling.ws3dinv.WSResponse(resp_fn=None, station_fn=None, wl_station_fn=None)[source]¶
class to deal with .resp file output by ws3dinv
Attributes
Description
n_z
number of vertical layers
period_list
list of periods inverted for
resp
- np.ndarray structured with keys:
station –> station name
- east –> relative eastern location in
grid
- north –> relative northern location in
grid
- z_resp –> impedance tensor array
of response with shape
(n_stations, n_freq, 4, dtype=complex)
*z_resp_err–> response impedance tensor error
resp_fn
full path to response file
station_east
location of stations in east direction
station_fn
full path to station file written by WSStation
station_names
names of stations
station_north
location of stations in north direction
units
[ ‘mv’ | ‘other’ ] units of impedance tensor
wl_sites_fn
full path to .sites file from Winglink
z_resp
impedance tensors of response with shape (n_stations, n_periods, 2, 2)
z_resp_err
impedance tensors errors of response with shape (n_stations, n_periods, 2, 2) (zeros)
Methods
Description
read_resp_file
read response file and fill attributes
Methods
read_resp_file
([resp_fn, wl_sites_fn, ...])read in data file
- class mtpy.modeling.ws3dinv.WSStartup(data_fn=None, initial_fn=None, **kwargs)[source]¶
read and write startup files
- Example
>>> import mtpy.modeling.ws3dinv as ws >>> dfn = r"/home/MT/ws3dinv/Inv1/WSDataFile.dat" >>> ifn = r"/home/MT/ws3dinv/Inv1/init3d" >>> sws = ws.WSStartup(data_fn=dfn, initial_fn=ifn)
Attributes
Description
apriori_fn
full path to a priori model file default is ‘default’
control_fn
full path to model index control file default is ‘default’
data_fn
full path to data file
error_tol
error tolerance level default is ‘default’
initial_fn
full path to initial model file
lagrange
starting lagrange multiplier default is ‘default’
max_iter
max number of iterations default is 10
model_ls
model length scale default is 5 0.3 0.3 0.3
output_stem
output file name stem default is ‘ws3dinv’
save_path
directory to save file to
startup_fn
full path to startup file
static_fn
full path to statics file default is ‘default’
target_rms
target rms default is 1.0
Methods
read_startup_file
([startup_fn])read startup file fills attributes
makes a startup file for WSINV3D.
- class mtpy.modeling.ws3dinv.WSStation(station_fn=None, **kwargs)[source]¶
read and write a station file where the locations are relative to the 3D mesh.
Attributes
Description
east
array of relative locations in east direction
elev
array of elevations for each station
names
array of station names
north
array of relative locations in north direction
station_fn
full path to station file
save_path
path to save file to
Methods
Description
read_station_file
reads in a station file
write_station_file
writes a station file
write_vtk_file
writes a vtk points file for station locations
Methods
from_wl_write_station_file
(sites_file, out_file)write a ws station file from the outputs of winglink
read_station_file
([station_fn])read in station file written by write_station_file
write_station_file
([east, north, ...])write a station file to go with the data file.
write_vtk_file
(save_path[, vtk_basename])write a vtk file to plot stations
- from_wl_write_station_file(sites_file, out_file, ncol=5)[source]¶
write a ws station file from the outputs of winglink
- mtpy.modeling.ws3dinv.cmap_discretize(cmap, N)[source]¶
Return a discrete colormap from the continuous colormap cmap.
cmap: colormap instance, eg. cm.jet. N: number of colors.
- Example
x = resize(arange(100), (5,100)) djet = cmap_discretize(cm.jet, 5) imshow(x, cmap=djet)
- mtpy.modeling.ws3dinv.computeMemoryUsage(nx, ny, nz, n_stations, n_zelements, n_period)[source]¶
compute the memory usage of a model
- mtpy.modeling.ws3dinv.estimate_skin_depth(res_model, grid_z, period, dscale=1000)[source]¶
estimate the skin depth from the resistivity model assuming that
delta_skin ~ 500 * sqrt(rho_a*T)