MTgeophysics / mtpy

Python toolbox for standard Magnetotelluric (MT) data analysis
GNU General Public License v3.0
147 stars 66 forks source link

abnormal mesh #164

Closed wilsonweijun closed 2 years ago

wilsonweijun commented 2 years ago

The Z-direction blocks are abnormal.

Expected Behavior

the mesh should be sensible. There are three option for padding methods. 1)extent1, 2)extent2),3)stretch, Coud you give more desriptions about these methods by words?

Current Behavior

the Model mesh:

MODEL FILE WRITTEN BY MTPY.MODELING.MODEM

67 149 72 0 LOGE 33710.000 10610.000 3340.000 1060.000 330.000 100.000 50.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.600 28.800 30.600 30.000 30.000 30.000 30.000 30.000 30.600 29.400 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.600 28.800 30.600 30.600 28.800 30.600 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 50.000 100.000 330.000 1060.000 3340.000 10610.000 33710.000 32820.000 10380.000 3280.000 1040.000 330.000 100.000 50.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 29.400 30.600 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 29.400 30.600 30.000 30.000 30.600 29.400 30.000 30.000 29.400 30.600 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.600 29.400 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.600 28.800 30.600 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 30.000 50.000 100.000 330.000 1040.000 3280.000 10380.000 32820.000 11.000 10.000 10.000 9.000 8.000 8.000 7.000 6.000 7.000 5.000 6.000 5.000 5.000 5.000 6.000 6.000 7.000 7.000 8.000 9.000 9.000 10.000 11.000 12.000 13.000 14.000 15.000 16.000 17.000 19.000 20.000 22.000 24.000 25.000 27.000 30.000 32.000 35.000 37.000 40.000 44.000 47.000 51.000 55.000 60.000 64.000 69.000 75.000 81.000 88.000 95.000 100.000 100.000 100.000 100.000 100.000 200.000 200.000 200.000 200.000 200.000 200.000 200.000 300.000 300.000 300.000 400.000 700.000 1000.000 1500.000 2300.000 3400.000

the thicknesses of Z-dericotn blocks are 11, 10, 9, 8,.... it is not descentding sequence . the thickness of the first earth layer is 5 m. the first layer thickness of air layers should be less than 5 meter, right?

Possible Solution

Steps to Reproduce (for bugs)

1. 2. 3. 4.


# Paste your code here
#
#
if not os.path.exists(workdir):
    os.mkdir(workdir)

edipath='D:/mtpywin/mtpy/HeiheEDIwithRot12deg'
edi_list=[os.path.join(edipath,ff) for ff in os.listdir(edipath) if(ff.endswith('.edi'))]

period_list=get_period_list(1.e-4,0.2,10)
do=Data(edi_list=edi_list,
        inv_mode='2',
        save_path=workdir,
        period_list=period_list,
        period_buffer=2,
        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
        model_epsg=4530)

do.write_data_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])

# Creat a model input file
from mtpy.modeling.modem import Model
workdir='D:/Manuscripts/Manuscript2020/ModemInv'
# create a Model object
mo=Model(station_locations=do.station_locations,
         cell_size_east=30,
         cell_size_north=30,
         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 cellls
         pad_num=6,         # number of cells outside station area before padding
         pad_stretch_v=1.5, # increasing factor in padding cells(vertical)
         pad_stretch_h=1.6, # increasing factor in padding cells(horizontal)
         n_air_layers=12,   # number of air layers
         res_initial_value=100,   # halfspace resistivity value for reference model (ohm-m)
         n_layers=60,       # total number of z layers,including air
         z1_layer=5,         # first layer thickness
         pad_method='extent1',   # method for calculating padding
         z_target_depth=4000   # approx. depth to bottom of core model
         )                     # (padding after this depth)
mo.make_z_mesh_new()   
mo.make_mesh()               
mo.write_model_file(save_path=workdir)

# add topography to the model and rewrite
topo_file='D:/DEM12.5m/heihe_AMT.asc'
mo.add_topography_to_model2(topo_file)
mo.write_model_file(save_path=workdir)
# project the stations on the topography
do.project_stations_on_topography(mo)

#Creat a covariance input file
from mtpy.modeling.modem import Covariance
co=Covariance()
co.smoothing_east=0.5
co.smoothing_north=0.5
co.smoothing_z=0.5
co.write_covariance_file(model_fn=mo.model_fn)

## Context
<!--- How has this issue affected you? What are you trying to accomplish? -->
<!--- Providing context helps us come up with a solution that is most useful in the real world -->
I hope to make a 3-D mesh to conduct AMT 3D inverison by using ModEM and MTPY.
## Your Environment
<!--- Include as many relevant details about the environment you experienced the bug in -->
  * Operating system: windows 10
  * MtPy version:  1.0 downloaded in 2020, with mesh_tools.py update from the patch 157.
  * Python version: 3.7.10
<!---if it is data visualization related, also provide-->
  * Matplotlib version:
  * Matplotlib backend (`print(matplotlib.get_backend())`):
<!---if it is graphical user interface (GUI) related-->
  * QT version:

**Installed Python Packages:**
use `pip freeze` or `conda list [-n ENVIRONMENT_NAME]` to list all the installed libraries.