|
| 1 | +""" |
| 2 | +Video Tutorial "code-along": Onlap relations |
| 3 | +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ |
| 4 | +""" |
| 5 | + |
| 6 | + |
| 7 | +# %% |
| 8 | +# This tutorial demonstrates step-by-step how to incorporate onlap relations to our geological models created with gempy. |
| 9 | +# It follows the Video tutorial series available on the [gempy YouTube channel](https://www.youtube.com/@GemPy3D). |
| 10 | +# Please follow the first and second part of the tutorials to learn the basics of modeling with gempy before diving into this tutorial. |
| 11 | + |
| 12 | +# %% |
| 13 | +# Video tutorial 11: Basic onlap scenario |
| 14 | +# """""""""""""""""""""""""""""""""""" |
| 15 | + |
| 16 | + |
| 17 | +# %% |
| 18 | +# .. raw:: html |
| 19 | +# |
| 20 | +# <iframe width="560" height="315" |
| 21 | +# src="https://www.youtube.com/embed/80QjnrFxubQ" |
| 22 | +# title="YouTube video player" |
| 23 | +# frameborder="0" |
| 24 | +# allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" |
| 25 | +# allowfullscreen> |
| 26 | +# </iframe> |
| 27 | +# |
| 28 | +# |
| 29 | + |
| 30 | +# %% |
| 31 | + |
| 32 | +# Required imports |
| 33 | +import gempy as gp |
| 34 | +import gempy_viewer as gpv |
| 35 | +import numpy as np |
| 36 | + |
| 37 | +# %% |
| 38 | + |
| 39 | +# Path to input data |
| 40 | +data_path = "https://raw.githubusercontent.com/cgre-aachen/gempy_data/master/" |
| 41 | +path_to_data = data_path + "/data/input_data/video_tutorials_v3/" |
| 42 | + |
| 43 | +# %% |
| 44 | + |
| 45 | +# Create instance of geomodel |
| 46 | +geo_model_onlap = gp.create_geomodel( |
| 47 | + project_name = 'tutorial_model_onlap_1', |
| 48 | + extent=[0,2000,0,1000,0,1000], |
| 49 | + resolution=[100,50,50], |
| 50 | + importer_helper=gp.data.ImporterHelper( |
| 51 | + path_to_orientations=path_to_data+"tutorial_model_onlap_1_orientations.csv", |
| 52 | + path_to_surface_points=path_to_data+"tutorial_model_onlap_1_surface_points.csv" |
| 53 | + ) |
| 54 | +) |
| 55 | + |
| 56 | + |
| 57 | +# %% |
| 58 | + |
| 59 | +# Map geological series to surfaces |
| 60 | +gp.map_stack_to_surfaces( |
| 61 | + gempy_model=geo_model_onlap, |
| 62 | + mapping_object={ |
| 63 | + "Young_Series": ("basin_fill_2", "basin_fill_1"), |
| 64 | + "Old_Series": ("basin_top", "basin_bottom") |
| 65 | + } |
| 66 | +) |
| 67 | + |
| 68 | +# Alternative way of mapping geological series to surfaces |
| 69 | +# gp.map_stack_to_surfaces( |
| 70 | +# gempy_model=geo_model_onlap, |
| 71 | +# mapping_object={ |
| 72 | +# "Young_Series": ("basin_fill_2", "basin_fill_1"), |
| 73 | +# "Onlap_Series": ("basin_top"), |
| 74 | +# "Old_Series": ("basin_bottom") |
| 75 | +# } |
| 76 | +# ) |
| 77 | + |
| 78 | +# %% |
| 79 | + |
| 80 | +# Display a basic cross section of input data |
| 81 | +gpv.plot_2d(geo_model_onlap, show_data=True) |
| 82 | + |
| 83 | +# %% |
| 84 | + |
| 85 | +# Compute a solution for the model |
| 86 | +gp.compute_model(geo_model_onlap) |
| 87 | + |
| 88 | +# %% |
| 89 | + |
| 90 | +# Display the result in 2d section |
| 91 | +gpv.plot_2d(geo_model_onlap, show_boundaries=False) |
| 92 | + |
| 93 | +# %% |
| 94 | + |
| 95 | +# Set the relation of the youngest group to Onlap |
| 96 | +from gempy_engine.core.data.stack_relation_type import StackRelationType |
| 97 | +geo_model_onlap.structural_frame.structural_groups[0].structural_relation = StackRelationType.ONLAP |
| 98 | + |
| 99 | +# %% |
| 100 | + |
| 101 | +# Display updated strucutral frame |
| 102 | +geo_model_onlap.structural_frame |
| 103 | + |
| 104 | +# %% |
| 105 | + |
| 106 | +# Compute a solution for the model |
| 107 | +gp.compute_model(geo_model_onlap) |
| 108 | + |
| 109 | +# %% |
| 110 | + |
| 111 | +# Display the result in 2d section |
| 112 | +gpv.plot_2d(geo_model_onlap, show_boundaries=False) |
| 113 | + |
| 114 | +# %% |
| 115 | +# Video tutorial 12: Advanced onlap - Subduction zone |
| 116 | +# """""""""""""""""""""""""""" |
| 117 | + |
| 118 | +# %% |
| 119 | +# .. raw:: html |
| 120 | +# |
| 121 | +# <iframe width="560" height="315" |
| 122 | +# src="https://www.youtube.com/embed/R-vUld4V-OQ" |
| 123 | +# title="YouTube video player" |
| 124 | +# frameborder="0" |
| 125 | +# allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" |
| 126 | +# allowfullscreen> |
| 127 | +# </iframe> |
| 128 | +# |
| 129 | +# |
| 130 | + |
| 131 | +# %% |
| 132 | + |
| 133 | +# Create instance of geomodel |
| 134 | +geo_model_subduction = gp.create_geomodel( |
| 135 | + project_name = 'tutorial_model_onlap_2', |
| 136 | + extent=[0,2000,0,1000,0,1000], |
| 137 | + resolution=[100,50,50], |
| 138 | + importer_helper=gp.data.ImporterHelper( |
| 139 | + path_to_orientations=path_to_data+"tutorial_model_onlap_2_orientations.csv?", |
| 140 | + path_to_surface_points=path_to_data+"tutorial_model_onlap_2_surface_points.csv?" |
| 141 | + ) |
| 142 | +) |
| 143 | + |
| 144 | +# %% |
| 145 | + |
| 146 | +# Display a basic cross section of input data |
| 147 | +gpv.plot_2d(geo_model_subduction) |
| 148 | + |
| 149 | +# %% |
| 150 | + |
| 151 | +# Map geological series to surfaces |
| 152 | +gp.map_stack_to_surfaces( |
| 153 | + gempy_model=geo_model_subduction, |
| 154 | + mapping_object={ |
| 155 | + "Top": ("continental_top"), |
| 156 | + "Continental_Series": ("continental_shallow", "continental_deep"), |
| 157 | + "Oceanic_Series": ("oceanic_top", "oceanic_bottom") |
| 158 | + } |
| 159 | +) |
| 160 | + |
| 161 | +# %% |
| 162 | + |
| 163 | +# Set the realtion of the youngest and second youngest group to Onlap |
| 164 | +geo_model_subduction.structural_frame.structural_groups[0].structural_relation = StackRelationType.ONLAP |
| 165 | +geo_model_subduction.structural_frame.structural_groups[1].structural_relation = StackRelationType.ONLAP |
| 166 | + |
| 167 | +# %% |
| 168 | + |
| 169 | +# Display updated structural frame |
| 170 | +geo_model_subduction.structural_frame |
| 171 | + |
| 172 | +# %% |
| 173 | + |
| 174 | +# Create a simple topography using numpy |
| 175 | + |
| 176 | +# Define grid spacing |
| 177 | +spacing = 20 |
| 178 | + |
| 179 | +# Generate grid |
| 180 | +x = np.arange(geo_model_subduction.grid.regular_grid.extent[0], geo_model_subduction.grid.regular_grid.extent[1] + spacing, spacing) |
| 181 | +y = np.arange(geo_model_subduction.grid.regular_grid.extent[2], geo_model_subduction.grid.regular_grid.extent[3] + spacing, spacing) |
| 182 | +X, Y = np.meshgrid(x, y) |
| 183 | + |
| 184 | +# Define elevation (z) based on x, creating a simple mountain range |
| 185 | +Z = np.ones_like(X) * 590 # Default elevation |
| 186 | +Z[(X >= 570) & (X < 1000)] = 590 + (200 * (X[(X >= 570) & (X < 1000)] - 600) / 400) |
| 187 | +Z[(X >= 1000) & (X < 1300)] = 810 - (250 * (X[(X >= 1000) & (X < 1300)] - 1000) / 300) |
| 188 | +Z[X >= 1300] = 540 |
| 189 | + |
| 190 | +# Flatten the data into (N,3) shape |
| 191 | +topography_points = np.vstack((X.ravel(), Y.ravel(), Z.ravel())).T |
| 192 | + |
| 193 | +# %% |
| 194 | + |
| 195 | +# Set topography from numpy array |
| 196 | +gp.set_topography_from_arrays(grid=geo_model_subduction.grid, xyz_vertices=topography_points) |
| 197 | + |
| 198 | +# %% |
| 199 | + |
| 200 | +# Compute a solution for the model |
| 201 | +gp.compute_model(geo_model_subduction) |
| 202 | + |
| 203 | +# %% |
| 204 | + |
| 205 | +# Display the result in 2d section |
| 206 | +gpv.plot_2d(geo_model_subduction, show_topography=True, show_boundaries=False) |
| 207 | + |
| 208 | +# %% |
| 209 | + |
| 210 | +# Display 3d plot of final model |
| 211 | +gpv.plot_3d(geo_model_subduction, show_topography=True, image=True) |
| 212 | + |
| 213 | +# sphinx_gallery_thumbnail_number = -1 |
| 214 | + |
0 commit comments