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.