/pyCSAMT

Python for Controlled Source Audio-frequency Magnetotellurics (CSAMT)

Primary LanguagePythonGNU Lesser General Public License v3.0LGPL-3.0

pyCSAMT : A Python open-source toolkit for Controlled Source Audio-frequency Magnetotellurics (CSAMT)

Documentation Status Build Status Requirements Status GitHub GitHub release (latest by date) DOI

Overview

  • Definition

    CSAMT is a geophysical method well-established as a resistivity exploration tool in deep geological structure detection. The method is broadly applied in diverse of exploration problems such as mineral , hydrocarbon, groundwater resources, as well as mapping the fault-zones etc.

  • Purpose

    The software contains bacics steps and improve CSAMT standard data processing and deals with OCCAM2D for modeling part. It also contains its inner database composed of geological structures and electrical properties of rocks, based on representative chart of Palacky (1988) and the rock and mineral property classification of Slichter and Telkes (1942) to generate a pseudo-stratigraphy log for drilling operations.

  • Note

    Actually pyCSAMT only works in far field and several outputs are provided for other external modeling softwares such as MTpy, OasisMontaj and GoldenSoftware.

Documentation

Licence

pyCSAMT is under GNU Lesser GPL version3 LGPLv3.

Units used

  • Frequency : [F] in Hz
  • Skin depth (sigma): sigma = 503 *sqrt([Rho]/[F]) in meters(m).
  • Apparent resistivy(Rho) : in Ω.m
  • E-field magnitude : [E]= microvolt/meter (muv/m)
  • H-field magnitude : [H] = gamma /amp
  • Impedance Tensor [Z] in 2*2 matrices : [Z] = [E]/[H]: km/s
  • Angle : Theta in degrees clockwise
  • Location coordinates ( X =N-S , Y = E-W) in m.
  • Coordinates in (UTM- Easting, Northing ) m.
  • Geomagnetic North : 0 degree azimuth
  • Step descent in m.
  • Input true resistivities in Ω.m

Available filters

  1. Trimming moving average (TMA) mostly used by Zonge International Engineering .
  2. Fixed-length-dipole moving average (FLMA) also used by Zonge International Engineering.
  3. Adaptative moving-average (AMA) based on idea of Torres-Verdin.
  4. MT Removal distorsion (dist) and static shift removal (ss) filters basically used to correct magnetotellurics (MT) data.

Plot inversion misfit and geostratigraphy misfit (misfit G)

To plot the misfit from measured data and the calculated inversion data, bring the occam response file (*.rep) and Occamlogfile (optional *.logfile) then run the script below:

  1. Plot some fitting curves of resistivity and phase inversion after applying on observed data the static shift correction.
>>> from pycsamt.modeling.occam2d import plotResponse 
>>> resPath =r'data/inversionFiles'                  # path to inversion files for each line
>>> _=plotResponse(data_fn =resPath,
...                 stations = ['S00', 'S04', 's08', 'S12'],  # sites to visualize 
...                  rms =['1.013', '1.451', '1.00', '1.069'], # rms of each line
...                  error_type ='resi' )

Click here to see the reference output.

  1. To plot the misfitof the model response from the FE algorithms:
>>> import os
>>> from pycsamt.modeling.occam2d import getMisfit 
>>> path_data ='data/occam2D'
>>> _=getMisfit(response_fn = os.path.join(path_data,'RESP17.resp'),
...         logfile=os.path.join(path_data, 'LogFile.logfile'), 
...          data_fn = path_data)

To see the output, click here.

  1. To evaluate the model errors misfit G between the the new resistivity model or stratigraphy models(NMs) from inversion models(CRMs), set misfit_G argument to True . Misfit G computation is the best way to see whether different layers with their corresponding resistivity values are misclassified or not. With few step of codes we can check the process:
>>> from pycsamt.geodrill.geocore import GeoStratigraphy
>>> inversion_files = {'model_fn':'data/occam2D/Occam2DModel', 
                       'mesh_fn': 'data/occam2D/Occam2DMesh',
                        "iter_fn":'data/occam2D/ITER17.iter',
                       'data_fn':'data/occam2D/OccamDataFile.dat'}
>>> resistivity_values =[10, 60, 70, 180, 1000,  3000, 7000]   # resistivity values of layers to map
>>> layer_names =['river water','sedimentary rocks', 'fracture zone',  'gravel', 'granite', 'igneous rocks','basement rocks' ]
>>> geosObj = GeoStratigraphy(**inversion_files,
...                      input_resistivities=resistivity_values, 
...                      input_layers=layer_names)
>>> geosObj.stratigraphyModel(kind='nm', misfit_G =False)           # 'nm':New Model

click here for reference output.

  • Note : For CSAMT data processing and some codes implementation, please refer to our wiki page.

Plot the pseudostratigraphic log

Once the geostratigraphic model is built, we just need to call the model and extract at each station its corresponding pseudostratigraphic log using the script below:

>>> from pycsamt.geodrill.geocore import GeoStratigraphy
>>> station ='S00'      # station to visualize 
>>> zoom =None          
>>> GeoStratigraphy.plotPseudostratigraphic(station =station, zoom = zoom )
 
~~~~~~~~~~~~~~~~~~~~~~~~~~~~[ PseudoStratigraphic Details: Station = S00 ]~~~~~~~~~~~~~~~~~~~~~~~~~~~~
------------------------------------------------------------------------------------------------------
|      Rank |            Stratum             |         Thick-range(m)         |     Thickness(m)     |
------------------------------------------------------------------------------------------------------
|        1. |         fracture zone          |         0.0 ----- 6.0          |         6.0          |
|        2. |             gravel             |         6.0 ----- 13.0         |         7.0          |
|        3. |            granite             |        13.0 ----- 29.0         |         16.0         |
|        4. |         igneous rocks          |        29.0 ----- 49.0         |         20.0         |
|        5. |         basement rocks         |        49.0 ----- 249.0        |        200.0         |
|        6. |         igneous rocks          |       249.0 ----- 289.0        |         40.0         |
|        7. |            granite             |       289.0 ----- 529.0        |        240.0         |
|        8. |         igneous rocks          |       529.0 ----- 699.0        |        170.0         |
|        9. |         basement rocks         |       699.0 ----- 999.0        |        300.0         |
------------------------------------------------------------------------------------------------------
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Survey Line: Occam2D files properties ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|model = Occam2DModel     |iter  = ITER17.iter      |mesh  = Occam2DMesh      |data  = OccamDataFile.dat|
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

click here for reference output.

It's possible to zoom the most interesting part of the log by setting the argument zoom, the running script above with zoom ='25%' gives the following ouput.

Credits

We use or link some third-party software (beside the usual tool stack: numpy, scipy, matplotlib) and are grateful for all the work made by the authors of these awesome open-source tools:

System requirements

  • Python 3.6+

Contributors

  1. Key Laboratory of Geoscience Big Data and Deep Resource of Zhejiang Province, School of Earth Sciences, Zhejiang University, China.

  2. Department of Geophysics, School of Geosciences and Info-physics, Central South University, China.

  3. Equipe de Recherche Géophysique Appliquée, Laboratoire de Geologie Ressources Minerales et Energetiques, UFR des Sciences de la Terre et des Ressources Minières, Université Félix Houphouët-Boigny, Cote d'Ivoire.