Make Structured 3D Mesh
A common input into 3D inversions is a structured mesh with equal cells within the station area and some padding cells. There are some tools in MTpy that can help design a mesh based on station locations, add topography or bathymetry, write to file specific for the inversion code.
[1]:
%matplotlib widget
[2]:
from pathlib import Path
from mtpy import MTData
from mtpy.modeling import StructuredGrid3D
1. Read Data File
In the previous notebook we created a data file, lets read that in and use the station locations to create a structured mesh for finite element codes like ModEM, WS3DINV, JIF3D.
[3]:
mtd = MTData()
mtd.from_modem(r"example_modem_data_file.dat")
2. Create a Model
Now that we have station locations we can create a 3D mesh.
[4]:
mesh = StructuredGrid3D()
mesh.station_locations = mtd.station_locations
mesh.center_point = mtd.center_point
print(mesh.center_point)
MT Location:
--------------------
Latitude (deg): 38.873786
Longitude (deg): -118.196245
Elevation (m): 0.0000
Datum crs: EPSG:4326
Easting (m): 396229.653
Northing (m): 4303450.561
UTM crs: EPSG:32611
Model Easting (m): 396229.653
Model Northing (m): 4303450.561
Model Elevation (m): 0.000
Profile Offset (m): 0.000
[ ]:
Horizontal Set Cell Sizes
We can set the lateral cell sizes in the east and north direction. There is a trade-off between cell size, computation speed, and resolution. Ideally you would like a few cells between each station to allow the inversion to find a better solution, but you don’t want too many that it will take ages to calculate.
The example data was collected at about 8 km station spacing, so lets try 2km cell sizes within the station area.
[5]:
mesh.cell_size_east = 2000
mesh.cell_size_north = 2000
Horizontal Padding cells
There are a few different padding cells.
pad_numis the number of padding cells of the same size as the station cells that extent outside the station area (default is 3)pad_eastandpad_northdescribe the number of padding cells between the edge of the station area and the desired extent, these are calculated on an exponential distance.
[6]:
mesh.pad_east = 10
mesh.pad_north = 10
[7]:
mesh.ew_ext = 250000
mesh.ns_ext = 250000
Vertical Cells
Because MT is a diffusive method we want to have a vertical grid that increases in size as a function of depth. There are a few ways to do this. Here we provide the first layer thickness z1_layer, the number of layers n_layers and the target depth z_target_depth. Then the cells will increase geometrically from the z1_layer to z_target_depth. These data are meant to image regional structures using broadband data, so the first layer should be thin and target depth around 50 km.
[8]:
mesh.z1_layer = 30
mesh.n_layers = 50
mesh.z_target_depth = 50000
Vertical Padding cells
We also need some padding cells in the vertical direction. We can set the geometric factor if needed.
[9]:
mesh.pad_z = 5
mesh.pad_stretch_v = 1.5
Make Mesh
Now we can create the mesh and see what it looks like.
[10]:
mesh.make_mesh()
Structured3DMesh Model Object:
--------------------
Number of stations = 59
Mesh Parameter:
cell_size_east: 2000
cell_size_north: 2000
pad_east: 10
pad_north: 10
pad_num: 3
z1_layer: 30
z_target_depth: 50000
n_layers: 50
n_air_layers: 0
res_initial_value: 100.0
Dimensions:
e-w: 61
n-s: 62
z: 51 (without 7 air layers)
Extensions:
e-w: 252000.0 (m)
n-s: 250000.0 (m)
0-z: 161492.0 (m)
--------------------
Note: in the plot below the station locations in the vertical profile show elevation, this is only if you want elevation included, it can be ignored for now. We will add topography later and project the stations properly.
[11]:
plot_mesh = mesh.plot_mesh(x_limits=(-mesh.ew_ext/2, mesh.ew_ext/2), y_limits=(-mesh.ns_ext/2, mesh.ns_ext/2))
26:05:08T22:35:19 | WARNING | line:49 |mtpy.modeling.plots.plot_mesh | _plot_topography | Cannot find topography information, skipping
26:05:08T22:35:19 | WARNING | line:89 |mtpy.modeling.plots.plot_mesh | _plot_topography_ax2 | Cannot find topography information, skipping
3. Add Topography
In some cases its important to model topography. You can download a DEM and convert it to an ESRI ASCII file or a geoTIFF (preferred). Or you can use the data elevations for a crude estimation of the elevation. You can get topography data at ETOPO. Be sure to download the icesurface layer.
You need to add some air cells to the model to include topograpy. If you do not then only bathymetry is projected onto the mesh, which is useful for areas near the coast.
airlayer_type identifies the method to add topography cells to the model
log_downplaces the top of the model at the top of the topography then increases cell sizes downwards and resets the mesh.log_upadded topography cells with increasing sizes upwards from sea level.constantadds equal thickness cells from the top of topography to mean elevation within the station area. (can be computationally expensive if there’s a lot of topography.
max_elev can be set to the maximum topography height in the model, which can be useful if there are large mountains between stations that don’t need to be modeled.
[12]:
mesh.n_air_layers = 40
[13]:
from mtpy_data import TOPO_EXAMPLE_FILE
mesh.add_topography_to_model(TOPO_EXAMPLE_FILE, airlayer_type="constant")
Read geotiff surface file
26:05:08T22:35:20 | INFO | line:1620 |mtpy.modeling.structured_mesh_3d | add_topography_to_model | Using constant air layer thickness of 30 m, which results in 45.0 air layers to cover the topography
[14]:
plot_topography = mesh.plot_mesh(fig_num=5, x_limits=(-mesh.ew_ext/4, mesh.ew_ext/4), y_limits=(-mesh.ns_ext/4, mesh.ns_ext/4), z_limits=(2000, -2000))
[15]:
mesh
[15]:
Structured3DMesh Model Object:
--------------------
Number of stations = 59
Mesh Parameter:
cell_size_east: 2000
cell_size_north: 2000
pad_east: 10
pad_north: 10
pad_num: 3
z1_layer: 30
z_target_depth: 50000
n_layers: 50
n_air_layers: 45
res_initial_value: 100.0
Dimensions:
e-w: 61
n-s: 62
z: 96 (without 7 air layers)
Extensions:
e-w: 252000.0 (m)
n-s: 250000.0 (m)
0-z: 162842.0 (m)
--------------------
5. Write to File
Now write to a file. This will be in ModEM style.
[16]:
mesh.to_modem(save_path=Path().cwd())
26:05:08T22:35:21 | INFO | line:830 |mtpy.modeling.structured_mesh_3d | to_modem | Wrote file to: c:\Users\peaco\OneDrive\Documents\GitHub\mtpy-v2\docs\source\notebooks\ModEM_Model_File.rho
6. Write Covariance File
Now that topography is in the model, we need to freeze air cells in the covariance matrix to avoid computational issues.
[17]:
from mtpy.modeling.modem import Covariance
[18]:
cov = Covariance()
cov.grid_dimensions = mesh.res_model.shape
cov.write_covariance_file(
save_path=Path().cwd(), res_model=mesh.res_model
)
26:05:08T22:35:21 | INFO | line:232 |mtpy.modeling.modem.convariance | write_covariance_file | Wrote covariance file to c:\Users\peaco\OneDrive\Documents\GitHub\mtpy-v2\docs\source\notebooks\covariance.cov
7. Project Stations onto Topography
We need to project the stations onto topography within the model.
Hint: With ModEM and WS3DINV it is also important that if topography is included then the stations should be at the center of the cell, otherwise you get some weird results.
[19]:
mtd.center_stations(mesh)
[20]:
mtd.project_stations_on_topography(mesh)
[24]:
mesh.station_locations = mtd.station_locations
mesh.plot_mesh(x_limits=(-mesh.ew_ext/4, mesh.ew_ext/4), y_limits=(-mesh.ns_ext/4, mesh.ns_ext/4), z_limits=(2000, -2000))
[24]:
Plotting PlotMesh
[21]:
mtd.to_modem(Path().cwd().joinpath("example_modem_data_with_topography.dat"))
26:05:08T22:35:27 | WARNING | line:539 |mtpy.modeling.modem.data | _check_for_too_small_errors | Found errors with values less than 0.02 in t_zx 1 times. Setting error as t_zx x 0.02.
26:05:08T22:35:27 | WARNING | line:539 |mtpy.modeling.modem.data | _check_for_too_small_errors | Found errors with values less than 0.02 in t_zy 1 times. Setting error as t_zy x 0.02.
26:05:08T22:35:28 | INFO | line:652 |mtpy.modeling.modem.data | write_data_file | Wrote ModEM data file to c:\Users\peaco\OneDrive\Documents\GitHub\mtpy-v2\docs\source\notebooks\example_modem_data_with_topography.dat
[21]:
ModEM Data Object:
Number of impedance stations: 59
Number of tipper stations: 59
Number of phase tensor stations: 59
Number of periods: 168
Period range (s):
Min: 0.0013021
Max: 2048
Rotation angle: 0.0
Data center:
Latitude: 38.8720 deg Northing: 4303248.1108 m
Longitude: -118.1925 deg Easting: 396551.9404 m
Datum epsg: 4326 UTM epsg: 32611
Elevation: 0.0 m
Impedance data: True
Tipper data: True
Inversion Mode: Full_Impedance, Full_Vertical_Components
[22]:
mtd.station_locations.head(5)
[22]:
| survey | station | latitude | longitude | elevation | datum_epsg | east | north | utm_epsg | model_east | model_north | model_elevation | profile_offset | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | data | gv100 | 38.611 | -118.535 | -1437.40 | 4326 | 0.0 | 0.0 | None | -29000.0 | -28000.0 | -1545.999 | 0.0 |
| 1 | data | gv101 | 38.595 | -118.351 | -1540.55 | 4326 | 0.0 | 0.0 | None | -13020.0 | -30020.0 | -1725.999 | 0.0 |
| 2 | data | gv102 | 38.594 | -118.277 | -1554.80 | 4326 | 0.0 | 0.0 | None | -7000.0 | -30020.0 | -1635.999 | 0.0 |
| 3 | data | gv103 | 38.585 | -118.202 | -1543.90 | 4326 | 0.0 | 0.0 | None | -1000.0 | -32020.0 | -1575.999 | 0.0 |
| 4 | data | gv104 | 38.596 | -118.137 | -1801.80 | 4326 | 0.0 | 0.0 | None | 5000.0 | -30020.0 | -1935.999 | 0.0 |
[23]:
mesh.station_locations = mtd.station_locations
p = mesh.plot_mesh(z_limits=(-2500, 2500))
[ ]: