Need to add the ability to read geotiff files to add topography to models.

Should use rasterio as that's already a dependency.

Should also think about other types of elevation files.

Trying to create PT & IV shapefiles from EDI directory

import os
from mtpy.utils.shapefiles import create_phase_tensor_shpfiles, create_tipper_shpfiles

edi_dir = "EDI"
save_dir_edi = "SHAPE"
if not os.path.exists(save_dir_edi):
    os.makedirs(save_dir_edi)
create_phase_tensor_shpfiles(edi_dir, save_dir_edi, proj="WGS84", ellipse_size=1000)
create_tipper_shpfiles(edi_dir, save_dir_edi)

Gives the following error

ImportError Traceback (most recent call last)
Cell In[12], line 2
1 import os
----> 2 from mtpy.utils.shapefiles import create_phase_tensor_shpfiles, create_tipper_shpfiles
4 # Define the path to your directory containing EDI files
5 edi_dir = "EDI"

File c:\Users\pmishra\Anaconda3\envs\mtpy\lib\site-packages\mtpy\utils\shapefiles.py:11
2 """
3 Create shape files for phase tensor ellipses.
4 https://pcjericks.github.io/py-gdalogr-cookbook/vector_layers.html#create-a-new-shapefile-and-add-data
(...)
8 @author: jrpeacock
9 """
10 import mtpy.modeling.modem
---> 11 from mtpy.utils.gis_tools import project_point_ll2utm, get_utm_zone
13 try:
14 from osgeo import ogr, gdal, osr

ImportError: cannot import name 'get_utm_zone' from 'mtpy.utils.gis_tools' (c:\Users\pmishra\Anaconda3\envs\mtpy\lib\site-packages\mtpy\utils\gis_tools.py)

Hi,
I am tiring to plot pseudo-section of phasetensor with tipper but I faced with this error:

C:\Users\Mohammad\.conda\envs\MTpy-v2\Lib\site-packages\scipy\stats\_stats_mstats_common.py:182: RuntimeWarning: invalid value encountered in scalar divide
  slope = ssxym / ssxm
C:\Users\Mohammad\.conda\envs\MTpy-v2\Lib\site-packages\scipy\stats\_stats_mstats_common.py:196: RuntimeWarning: invalid value encountered in sqrt
  t = r * np.sqrt(df / ((1.0 - r + TINY)*(1.0 + r + TINY)))
C:\Users\Mohammad\.conda\envs\MTpy-v2\Lib\site-packages\scipy\stats\_stats_mstats_common.py:199: RuntimeWarning: invalid value encountered in scalar divide
  slope_stderr = np.sqrt((1 - r**2) * ssym / ssxm / df)
Traceback (most recent call last):
  File "C:\Users\Mohammad\Desktop\MTpy_V2\2.py", line 11, in <module>
    plot_pt_pseudosection = dgh_mt_data.plot_phase_tensor_pseudosection(
                            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\Mohammad\.conda\envs\MTpy-v2\Lib\site-packages\mtpy\core\mt_data.py", line 1116, in plot_phase_tensor_pseudosection
    return PlotPhaseTensorPseudoSection(mt_data=self, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\Mohammad\.conda\envs\MTpy-v2\Lib\site-packages\mtpy\imaging\plot_phase_tensor_pseudosection.py", line 101, in __init__
    self.plot()
  File "C:\Users\Mohammad\.conda\envs\MTpy-v2\Lib\site-packages\mtpy\imaging\plot_phase_tensor_pseudosection.py", line 340, in plot
    self._get_profile_line()
  File "C:\Users\Mohammad\.conda\envs\MTpy-v2\Lib\site-packages\mtpy\imaging\mtplot_tools\base.py", line 655, in _get_profile_line
    mt_obj.project_onto_profile_line(
  File "C:\Users\Mohammad\.conda\envs\MTpy-v2\Lib\site-packages\mtpy\core\mt_location.py", line 452, in project_onto_profile_line
    raise ValueError(
ValueError: utm_crs is None, cannot project onto profile line.

my script is:

from mtpy import MTCollection
mc = MTCollection()
mc.open_collection(r'C:\Users\Mohammad\Desktop\work\nasrabad_data\nasrabad_data\15\15.h5')
#plot_strike_all = mc.plot_strike()
mc.utm_crs = 32640
mc.working_dataframe = mc.master_dataframe[mc.master_dataframe.station.str.startswith("15_010-0000")]
dgh_mt_data = mc.to_mt_data()
plot_pt_pseudosection = dgh_mt_data.plot_phase_tensor_pseudosection( 
    plot_tipper="yri", 
    profile_reverse=False,
    y_limits=(10**2, 1)
)

Need to add more options when converting model files to other formats for viewing.

Options needed:

• xy coordinates [ lat/lon | easting/northing] and need to be able to specify units and coordinate system [ utm | datum ]
• z coordinates [ positive up | positive down ]
• need to able to specify a clip for all 3 dimensions.

The induction vectors are currently plotted in german convention where reals are pointing away from a good conductor, the default should be parkinson convention to point towards a good conductor. Just need to set arrow.arrow_direction = -1 as the default for plotting.

Expected Behavior

If an edi file does not have impedance, but only tipper, for a station that does not have electric dipoles, the plotter functions crash.

Current Behavior

from mtpy import MT, MTData
md = MTData()
md.add_station(list_of_paths_to_edi_files, survey="survey_name")
p = md.plot_mt_response("survey_name.001")
p = md.plot_phase_tensor_map()

Bascially, it looks like impedance values are being looked for but not found.

Maybe instead of letting Z be None it could be an empty iterable?

Need to fix converting to geographic coordinates when using StructuredGrid3D.to_vtk() and MTData.to_vtk.

data needs a coordinate system set ne+, en-

model needs to be tested cause it doesn't come out correctly.

Thanks for the good work you are doing. I keep getting below error when outputting convariace_topo.cov in jupyter.

TypeError Traceback (most recent call last)
Cell In[40], line 2
1 topo_covariance = Covariance(grid_dimensions=model_mesh.res_model.shape)
----> 2 topo_covariance.write_covariance_file(
3 Path().cwd().joinpath("modeling","musgraves_covariance_topo.cov"),
4 model_fn=model_mesh.model_fn,
5 )

TypeError: write_covariance_file() got an unexpected keyword argument 'model_fn'

Update documentation on readthedocs.

Should add ability to download transfer functions from various clients that host transfer functions.

The mtpy-v2 doesn't seem to have Occam2D Plot_Model functionality (for .iter file) like version 1. Also, I couldn't find any new function for plotting the resistivity model. I tried to use the following as per mtpy-v2 documentation:

import mtpy.modeling.occam2d as occam2d
itfn = r"C:/Users/kaush/Desktop/LD_sensitibility/ITER15.iter"
ocm = occam2d.Model(itfn)
ocm.read_iter_file()

But I get the error: module 'mtpy.modeling.occam2d' has no attribute 'Model'
When I change the Model to model (using lowercase m instead, the error is : 'module' object is not callable.
Also, I am not sure whether it will plot the resistivity model.
Please let me know of any solution.

Regards,
Kaushik Sarker

I am trying to add edis from my project (in Finland) to MTCollection

from pathlib import Path
from mtpy import MTCollection
tf_path=Path(r"../path_to_my_edis")
tf_path.exists()
with MTCollection() as mc:
    mc.open_collection(tf_path.joinpath("MyProject_tf.h5"))
    mc.add_tf(mc.make_file_list(tf_path)) 

Error

24:02:20T18:09:26 | ERROR | line:203 |mt_metadata.transfer_functions.io.tools | get_nm_elev | Input values (latitude=63.74171066666667, longitude=26.25233077777778) could not be found on US National Map.

Is there a way to read elevation from EDIs or perhaps define them in an external file?

when to create data for 3D inversion I faced with error:
interp_mt_data.calculate_rel_locations() AttributeError: 'MTData' object has no attribute 'calculate_rel_locations'

Should add ability to make NetCDF files from model files.

Looks like GIS support could be from rioxarray and xarray for NetCDF files.

Hello, I am running the following code in mtpy_v2 to generate ModEM inversion data.

-- coding: utf-8 --

"""
Created on Fri Nov 26 15:11:11 2021

@author: kaush
"""

import os
import mtpy.modeling.modem as modem
from mtpy.modeling.modem import Data
#from mtpy.core.edi import Edi
#from mtpy.utils.calculator import get_period_list
import numpy as np

path to save to

workdir = r'C:/tmp/ModEM_inv/test/test3D'

make the save path if it doesn't exist

if not os.path.exists(workdir):
os.mkdir(workdir)

path containing edi files

edipath = r'C:/Users/kaush/Desktop/mtpy-v1.0/mtpy/examples/data/test_3D'
#the number of files are the number of stations. Each file=each station

find all files in edipath

edi_list = [os.path.join(edipath,ff) for ff in os.listdir(edipath) if
(ff.endswith('.edi'))]

define period list to interpolate data onto

MTPy will not extrapolate outside the data period range

this code gets 4 periods per decade

#period_list = get_period_list(1e-3,1e3,1)
period_list = np.array([3.174e-2, 3.662e-2, 4.395e-2, 5.371e-2,
6.348e-2, 7.324e-2, 8.789e-2, 1.074e-1, 1.270e-1, 1.465e-1, 1.758e-1,
2.148e-1, 2.539e-1, 2.930e-1, 3.516e-1, 4.297e-1, 5.078e-1, 5.859e-1,
7.031e-1, 8.594e-1, 1.016e+0, 1.172e+0, 1.406e+0, 1.719e+0, 2.031e+0,
2.344e+0, 2.813e+0, 3.438e+0])

create a data object

do = Data(edi_list=edi_list,
inv_mode = '2', # invert for Z + tipper ('2' means Z only)
save_path=workdir,
period_list=period_list,
period_buffer=2., # buffer factor to determine how far to interpolate
# between points. 2 means interpolate by a factor of
# 2. e.g. if the data is at 10s the code will
# interpolate only from 5 to 20 s.

# #debugging error code 
#       error_type_z=np.array([['floor_egbert','floor_egbert'], ['floor_egbert','floor_egbert']]), # add floor to apply it as an error floor
                                                            
#       error_value_z=np.array([[1,4.5], [4.5,1]]), # can supply a 2 x 2 array for each component or a single value
# #debugging error code end

#original code error
      error_type_z='floor_egbert', # error type (egbert is % of sqrt(zxy*zyx))
      error_value_z=10., # error floor in percent
      error_type_tipper = 'floor_abs', # type of error to set in tipper,

floor_abs is an absolute value set

as a floor

      error_value_tipper=0.001,
      #original error code ending
      station_locations=modem.Stations().get_station_locations(edi_list),

     # station_locations=modem.station().get_station_locations(edi_list)


      model_epsg=2378
     
      #model_epsg=26996#2274 #NMSZ #28354 # model epsg, currently set to utm zone 54.

See http://spatialreference.org/ to

find the epsg code for your projection

    )

write data file, setting elevation to zero as we haven't added topography

do.write_data_file()

save epsg code and centre position to a text file

np.savetxt(os.path.join(workdir,'center_position.txt'),
np.array([do.center_point['east'],
do.center_point['north']]))
#np.savetxt(os.path.join(workdir,'epsg.txt'),[do.model_epsg])

############ Model ##################################
from mtpy.modeling.modem import Model

create a Model object

mo = Model(station_locations=do.station_locations,
cell_size_east=8000,
cell_size_north=8000,
pad_north=7, # number of padding cells N and S
pad_east=7,# number of padding cells E and W
pad_z=6, # number of vertical padding cells
pad_num=6, # number of cells outside station area before padding
pad_stretch_v=1, # increase factor in padding cells (vertical)
pad_stretch_h=1, # increase factor in padding cells (horizontal)
n_air_layers = 1, # number of air layers
res_initial_value=100, # halfspace resistivity value
# for reference model (ohm-m)
n_layers=10, # total number of z layers, including air
z1_layer=1000, # first layer thickness
pad_method='stretch', # method for calculating padding
z_target_depth=15000 # approx. depth to bottom of core model
# (padding after this depth)
)
mo.make_mesh()
mo.write_model_file(save_path=workdir)

from mtpy.modeling.modem import Covariance

create a covariance file

co = Covariance()
co.smoothing_east = 0.6
co.smoothing_north = 0.6
co.smoothing_z = 0.6
co.write_covariance_file(model_fn=mo.model_fn)

I am getting the following message: Invalid projection: EPSG:None: (Internal Proj Error: proj_create: crs not found). I am using data from New Madrid Seismic Zone. I can't find any specific epsg number or utm zone for the region. Please let me know about the problem.

Regards,
Kaushik Sarker

Need to add support for outputting phase tensor information to modeling data files.

reading of modem data file causes index error

md.read_data_file(fname)
File ~\Anaconda3\envs\mtpy-v2\lib\site-packages\mtpy\modeling\modem\data.py:1094 in read_data_file
n_periods, n_stations = self._read_header(
File ~\Anaconda3\envs\mtpy-v2\lib\site-packages\mtpy\modeling\modem\data.py:891 in _read_header
self._parse_header_line(head_line, inv_mode)
File ~\Anaconda3\envs\mtpy-v2\lib\site-packages\mtpy\modeling\modem\data.py:930 in _parse_header_line
value = item_list[1].replace("%", "").split()[0]
IndexError: list index out of range

Current Behavior

The first line of all of our ModeEM data files looks like this:
# ModEM impedance responses for Description:
item_list will contain ['description', ''] for each of our files and than fail

Your Environment

• MtPy version: 2.0.5

I have a MT collection that has no survey id and to downsample frequencies, I needed to add the survey id because the keys for MT collection are survey_id.station_id and having '0' for survey id meant no string and problem. I got this error:

I have used this code to downsample and add the id:

# with MTCollection() as mc:
     mc.open_collection(r"C:\Users\aivazpourporgous\Documents\MT\CCA\Inversion\CCA_0\edi\lin0_resampled_frequency.h5")
     for filename in edifiles:
        mt_object = MT()
        mt_object.read(filename)
        #mt_object.survey = "downsampled"
        new_frequency = np.logspace(-3, 3, 20)
        new_mt_obj = mt_object.interpolate(new_frequency, f_type="frequency")
        new_mt_obj.survey = "downsampled"
        mc.add_tf(new_mt_obj)
#
#

Now I have two entities for each station one unknown survey id and one downsampled. I added the id to the original mt_object before making new_mt_obj but result is same.

I'm trying to create a MTCollection for a large data set (> 1000 edi files). Using a loop with mt_object.read(fn) / md.add_station(mt_object), it takes > 5 hours. It slows down with iteration, as if there is some reinitialization/copy every time the add function is called.

Is this a feature or is there perhaps a more efficient bulk read function?

Add ability to output GMT files for plotting phase tensors, induction arrows, model maps.

Copying an issue posted in the old mtpy repository that I am also having. I want to create my own MT_collection with my MT_data. When I look for how to create it, the link created for the (mtpy_data) at below, does not open.

https://github.com/MTgeophysics/mtpy_data

Test data can be found at mtpy-data. Installation instructions are included in README.

from mtpy_data import GRID_LIST, PROFILE_LIST

Initiate MTCollection
from pathlib import Path
from mtpy import MT, MTCollection

2023-04-14 11:48:00,318 [line 135] mth5.setup_logger - INFO: Logging file can be found C:\Pathway\debug.log

mtc = MTCollection()
mtc.open_collection(Path().cwd().joinpath("test_mt_collection.h5"))