"""
Description:
With an input edi_file_folder and a list of period index,
generate a profile using occam2d module,
then plot the Penetration Depth profile at the given periods vs the stations locations.
Usage:
python mtpy/imaging/penetration_depth2d.py /path2/edi_files_dir/ period_index_list
python mtpy/imaging/penetration_depth2d.py.py examples/data/edi2/ 0 1 10 20 30 40
Author: fei.zhang@ga.gov.au
Date: 2017-01-23
"""
import glob
import os
import sys
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import click
from mtpy.imaging.penetration import get_penetration_depth, load_edi_files, Depth2D
# mpl.rcParams['lines.linewidth'] = 2
# mpl.rcParams['lines.color'] = 'r'
# mpl.rcParams['figure.figsize'] = [20, 10]
import mtpy.core.mt as mt
import mtpy.modeling.occam2d_rewrite as occam2d_new
from mtpy.utils.mtpylog import MtPyLog
# get a logger object for this module, using the utility class MtPyLog to
# config the logger
_logger = MtPyLog.get_mtpy_logger(__name__)
# logger =
# MtPyLog(path2configfile='logging.yml').get_mtpy_logger(__name__) #
# specific
# use the Zcompotent=[det, zxy, zyx]
def plot2Dprofile(edi_dir, period_index_list=None, zcomponent='det',
edi_list=None, tick_params={}, save=False,savepath=None,**kwargs):
#edi_dir = "/Softlab/Githubz/mtpy2/tests/data/edifiles/"
# edi_dir="E:/Githubz/mtpy2/tests/data/edifiles/"
# edi_dir=r"E:\Githubz\mtpy2\examples\data/edi2"
# 1 get a list of edi files, which are suppose to be in a profile.
edifiles = glob.glob(os.path.join(edi_dir, '*.edi'))
_logger.debug("edi files: %s", edifiles)
edis = load_edi_files(edi_dir, file_list=edi_list)
plot = Depth2D(edis, period_index_list, zcomponent)
plot.plot(tick_params, **kwargs)
if save:
if os.path.isdir(savepath):
savepath == os.path.join(savepath,'Depth2D.png')
if savepath is not None:
plot._fig.savefig(savepath)
else:
savepath = os.path.join(edi_dir,'Depth2D.png')
plot.show()
[docs]def barplot_multi_station_penentration_depth(
edifiles_dir, per_index=0, zcomponent='det'):
"""
A simple bar chart plot of the penetration depth across multiple edi files (stations),
at the given (frequency) per_index. No profile-projection is done in this funciton.
:param edifiles_dir: a list of edi files, or a dir of edi
:param per_index: an integer smaller than the number of MT frequencies in the edi files.
:return:
"""
if os.path.isdir(edifiles_dir):
edi_dir = edifiles_dir # "E:/Githubz/mtpy2/tests/data/edifiles/"
edifiles_dir = glob.glob(os.path.join(edi_dir, '*.edi'))
_logger.debug(edifiles_dir)
else:
# Assume edifiles_dir is [a list of edi files]
pass
scale_param = np.sqrt(1.0 / (2.0 * np.pi * 4 * np.pi * 10 ** (-7)))
# per_index=0,1,2,....
periods = []
depths = []
stations = []
mt_obj_list = [mt.MT(afile) for afile in edifiles_dir]
(stations, periods, depths, _) = get_penetration_depth(
mt_obj_list, int(per_index), whichrho=zcomponent)
# the attribute Z
# zeta = mt_obj.Z
#
# if per_index >= len(zeta.freq):
# raise Exception("Error: input period index must be less than number of freqs in zeta.freq=%s",len(zeta.freq))
#
# per = 1.0 / zeta.freq[per_index]
# periods.append(per)
# penetration_depth = -scale_param * np.sqrt(zeta.resistivity[per_index, 0, 1] * per)
#
# depths.append(penetration_depth)
# stations.append(mt_obj.station)
#plt.plot(app_resis, color='b', marker='o')
index = np.arange(len(depths))
plt.bar(index, depths, color='#000000')
# plt.xaxis.tick_top()
# plt.set_xlabel('X LABEL')
# plt.xaxis.set_label_position('top')
plt.xlabel(
'Penetration Depth Across Stations, for MT period= %6.5f Seconds' %
periods[0], fontsize=16)
plt.ylabel('Penetration Depth (m)', fontsize=16)
# plt.title('Penetration Depth profile for T=??')
bar_width = 0.4
plt.xticks(
index + bar_width / 2,
stations,
rotation='horizontal',
fontsize=14)
plt.legend()
# plt.tight_layout()
plt.gca().xaxis.tick_top()
plt.show()
# Check that the periods are the same value for all stations
return (stations, depths, periods)
# =============================================================================================
# Example Usage:
# python mtpy/imaging/penetration_depth2d.py examples/data/edi_files/ 1 10 20 30
# python mtpy/imaging/penetration_depth2d.py tests/data/edifiles/ 0 1 10 20 30 40 50 59
# python mtpy/imaging/penetration_depth2d.py examples/data/edi2/ 0 1 10 20 30 40
# =============================================================================================
if __name__ == "__main__old":
if len(sys.argv) < 2:
print(("Usage: %s edi_dir" % sys.argv[0]))
print ("python examples/penetration_depth2d.py tests/data/edifiles/ 0 1 10 20 30 40 50 59")
sys.exit(1)
elif os.path.isdir(sys.argv[1]):
edi_dir = sys.argv[1] # the first argument is path2_edi_dir
# the second,.... will be period index list
period_index_list = sys.argv[2:]
print(("period_index_list = {}".format(period_index_list)))
# the rho zcomponent can be det, zxy zyx
plot2Dprofile(edi_dir, period_index_list, zcomponent='det')
perindex = int(sys.argv[2])
# barplot_multi_station_penentration_depth(edi_dir, per_index=perindex)
# #, zcomponent='zxy')
else:
print("Please provide an edi directory and period_index_list")
# =============================================================================================
# Command line wrapper
# =============================================================================================
@click.command()
@click.option('-i','--input',type=str,default='examples/data/edi_files',help='directory or edsi data files')
@click.option('-p','--period_list',type=str,default="0 1 10 20 30 40",help='Periods seperated by space')
def plot_penetration_image(input,period_list):
if os.path.isdir(input):
period_index_list = period_list.split(' ')
plot2Dprofile(input, period_index_list, zcomponent='det')
else:
print("Please provide an edi directory !")
if __name__ == '__main__':
plot_penetration_image()