Modeling¶
Module ModEM¶
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
- rp_tuple : np.ndarray (freq, res, res_err, phase, phase_err)
- with res, phase having shape (num_freq, 2, 2).
- edi_file : string
- full path to edi file to be modeled.
- save_path : string
- path to save the file, if None set to dirname of station if edipath = None. Otherwise set to dirname of edipath.
- thetar : float
- 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_err : float
errorbar for resistivity values. Can be set to ( default = ‘data’):
- ‘data’ for errorbars from the data
- percent number ex. 10 for ten percent
- phase_err : float
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
()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
()redraw plot if parameters were changed save_figure
(save_fn[, file_format, …])save_plot will save the figure to save_fn. update_plot
()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 -
write_startup_file
(save_path=None, **kwargs)[source]¶ Make a 1D input file for Occam 1D
- savepath : full path to save input file to, if just path then
- saved as savepath/input
- model_fn : full path to model file, if None then assumed to be in
- savepath/model.mod
- data_fn : full 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_rho : starting resistivity value on linear scale.
- default = 100
description : description of the inversion.
- start_lagrange : starting Lagrange multiplier for smoothness.
- default = 5
start_rough : starting roughness value. default = 1E7
- debuglevel : something to do with how Fortran debuggs the code
- Almost always leave at default = 1
- start_iter : the 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 - data : is 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_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 * data * data_list * freq * station_list * station_locations 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 * data * data_list * freq * station_list * station_locations 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_mesh
()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_model
()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
()redraw plot if parameters were changed save_figure
(save_fn[, file_format, …])save_plot will save the figure to save_fn. update_plot
()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
get_misfit
()compute misfit of MT response found from the model and the data. plot
()plot pseudo section of data and response if given redraw_plot
()redraw plot if parameters were changed save_figure
(save_fn[, file_format, …])save_plot will save the figure to save_fn. update_plot
()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
()redraw plot if parameters were changed save_figure
(save_fn[, file_format, …])save_plot will save the figure to save_fn. update_plot
()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
()redraw plot if parameters were changed save_figure
(save_fn[, file_format, …])save_plot will save the figure to save_fn. update_plot
()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
()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_profile
()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. build_regularization
()Builds larger boxes around existing mesh blocks for the regularization. get_num_free_params
()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: * binding_offset * mesh_fn * model_columns * model_rows * prejudice_fn * statics_fn 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
get_misfit
()compute misfit of MT response found from the model and the data. plot
()plot pseudo section of data and response if given redraw_plot
()redraw plot if parameters were changed save_figure
(save_fn[, file_format, …])save_plot will save the figure to save_fn. update_plot
()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
()redraw plot if parameters were changed save_figure
(save_fn[, file_format, …])save_plot will save the figure to save_fn. update_plot
()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
()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_files
()read in the files to get appropriate information redraw_plot
()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
()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 residualecolor = ‘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 scalexpad = 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
()redraw plot if parameters were changed save_figure
(save_fn[, file_format, …])save_plot will save the figure to save_fn. update_plot
()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_station_grid_locations
()get the grid line on which a station resides for plotting on_key_press
(event)on a key press change the slices plot
()plot: east vs. read_files
()read in the files to get appropriate information redraw_plot
()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
build_data
()Builds the data from .edi files to be written into a data file compute_errors
()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_model_to_int
()convert the resistivity model that is in ohm-m to integer values corresponding to res_list make_mesh
()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_model_file
()read in a model file as x-north, y-east, z-positive down write_vtk_file
(save_fn)
-
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_model_to_int
()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: -a radio dial for depth slice -radio dial for resistivity value read_file
()reads in initial file or model file and set attributes: -resmodel -northrid -eastrid -zgrid -res_list if initial file rect_onselect
(eclick, erelease)on selecting a rectangle change the colors to the resistivity values redraw_plot
()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
(label)-
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 write_startup_file
()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)