Can't make ModEM file due to not knowing epsg.

[kaushik67]

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

[oaazeved]

Hi there! Please edit your comment to make it more readable; you can surround a block of code with the characters ``` on either side to prevent your code comments and other characters from being parsed as Markdown by Github.

Regarding an appropriate EPSG number, consider referring to epsg.io; a quick glance suggests to me that EPSG:26996 may not be a bad choice.

[kujaku11]

@kaushik67 I would have a look at https://github.com/kujaku11/mt_examples/blob/main/notebooks/mtpy/14_modem_outputs.ipynb, it looks like you might be running an older script for v1.