{ "cells": [ { "cell_type": "markdown", "id": "worldwide-contact", "metadata": {}, "source": [ "\n", "\n", "#
From Maps to Models - Tutorials for structural geological modeling using GemPy and GemGIS
\n", "\n", "# Model 2 - Folded Layers\n", "\n", "The second notebook in this tutorial series builds upon the [first notebook](model1_Horizontal_Layers.ipynb) ([notebook on Github](https://nbviewer.org/github/cgre-aachen/gemgis_data/blob/main/notebooks/01_basic_modeling/model1_Horizontal_Layers.ipynb)) where parallel horizontal layers were modeled. This notebook illustrates now how to create a simple model of folded layers in `GemPy`. The model consists of four parallel but folded layers plus basement layers and has an extent of 1000 m by 1000 m with a vertical extent of 600 m. No faulted or truncated layers are present in the model. \n", "\n", "If you have not gone through the introduction notebook for the course, please check it out: [Introduction Notebook](../00_introduction_to_structural_modeling.ipynb) ([notebook on Github](https://nbviewer.org/github/cgre-aachen/gemgis_data/blob/main/notebooks/00_introduction_to_structural_modeling.ipynb))\n", "\n", "\n", "
\n", "In this tutorial, you will learn the following:
\n", "- Get a better undestanding of what orientations mean in GemPy and how to plot them using mplstereonet
\n", "- How to build a simple model consisting of folded layers belonging to one Series
\n", "\n", "
\n", "\n", "## Contents\n", "\n", "1. [Installing GemPy](#installing-gempy)\n", "2. [Importing Libraries](#importing-libraries)\n", "3. [Data Preparation](#data-preparation)\n", " 1. [Importing Interface Points](#importing-interface-points)\n", " 2. [Importing Orientations](#importing-orientations)\n", " 1. [Visualizing the Orientations using mplstereonet](#visualizing-the-orientations-using-mplstereonet)\n", "4. [GemPy Model Calculation](#gempy-model-calculation)\n", " 1. [Creating the GemPy Model](#creating-the-gempy-model)\n", " 2. [Data Initiation](#data-initiation)\n", " 3. [Inspecting the Surfaces](#inspecting-the-surfaces)\n", " 4. [Inspecting the Input Data](#inspecting-the-input-data)\n", " 5. [Map Stack to Surfaces](#map-stack-to-surfaces)\n", " 6. [Plotting Input Data in 2D](#plotting-the-input-data-in-2d)\n", " 7. [Plotting Input Data in 3D](#plotting-the-input-data-in-3d)\n", " 8. [Setting the Interpolator](#setting-the-interpolator)\n", " 9. [Computing the Model](#computing-the-model)\n", "5. [Model Visualization and Post-Processing](#model-visualization-and-post-processing)\n", " 1. [Visualizing Cross Sections of the Computed Model](#visualizing-cross-sections-of-the-computed-model)\n", " 2. [Visualizing the computed model in 3D](#visualizing-the-computed-model-in-3d)\n", "6. [Conclusions](#conclusions)\n", "7. [Outlook](#outlook)\n", "8. [Licensing](#licensing)\n", "\n", "\n", "\n", "" ] }, { "cell_type": "markdown", "id": "1df9cae6", "metadata": {}, "source": [ "The input data is provided as already prepared CSV-files and will be loaded as Pandas `DataFrames`. The orientation file contains the orientation of the layers with with the locations of the orientations (`X`, `Y`, `Z`), the associated `formation` plus three more columns, `dip`, `azimuth`, and `polarity` indicating the dip the layer and the dip direction (azimuth, see Figure below) in degrees. The dip varies from 0° for horizontal layers to 90° for vertical layers. The azimuth varies from 0° (N) via 180° (S) to 360° (N). Here, we only provide the orientations for one layer. This will be explained later on in more detail. The `polarity` value is mostly set to 1.\n", "\n", "\n", "\n", "By CrunchyRocks, after Karla Panchuck - https://openpress.usask.ca/physicalgeology/chapter/13-5-measuring-geological-structures/, CC BY 4.0, https://commons.wikimedia.org/w/index.php?curid=113554289\n", "\n", "\n", "Orientations are, apart from interface points, the second important type of input data for `GemPy`. By providing orientations, you will provide a constraint for `GemPy` for the gradient of the scalar field. This is illustrated in the figure below where the colored lines represent the calculated scalar field, the arrows represent orientations, where the azimuth of the orientation defines the direction of the greatest gradient of the scalar field and the dip represents the \"magnitude\" of the gradient (the higher the dip, the higher the gradient, the narrower the isolines). The interfaces points (blue and red dots) represent two parallel folded layers that are located on an isoline of the scalar field. The isoline of the scalar field thus represents the boundary between two layers and would be termed a surface.\n", "\n", "\n", "Source: de la Varga et al. (2019)" ] }, { "cell_type": "markdown", "id": "solid-hormone", "metadata": {}, "source": [ "\n", "\n", "# Installing GemPy\n", "\n", "If you have not installed `GemPy` yet, please follow the [installation instructions](https://docs.gempy.org/installation.html). If you encounter any issues, feel free to open a new discussion at [GemPy Discussions](https://github.com/cgre-aachen/gempy/discussions). If you encounter an error in the installation process, feel free to also open an issue at [GemPy Issues](https://github.com/cgre-aachen/gempy/issues). There, the `GemPy` development team will help you out. " ] }, { "cell_type": "markdown", "id": "artificial-customs", "metadata": {}, "source": [ "\n", "\n", "# Importing Libraries\n", "\n", "For this notebook, we need the `pandas` library for the data preparation, `matplotlib` for plotting, the `mplstereonet` library for plotting orientation data and of course the `gempy` library. Any warnings that may appear can be ignored for now. " ] }, { "cell_type": "code", "execution_count": 1, "id": "2f498618", "metadata": { "ExecuteTime": { "end_time": "2022-04-02T12:33:29.872541Z", "start_time": "2022-04-02T12:33:27.083564Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "WARNING (theano.configdefaults): g++ not available, if using conda: `conda install m2w64-toolchain`\n", "C:\\Users\\ale93371\\Anaconda3\\envs\\gempy_new8\\lib\\site-packages\\theano\\configdefaults.py:560: UserWarning: DeprecationWarning: there is no c++ compiler.This is deprecated and with Theano 0.11 a c++ compiler will be mandatory\n", " warnings.warn(\"DeprecationWarning: there is no c++ compiler.\"\n", "WARNING (theano.configdefaults): g++ not detected ! Theano will be unable to execute optimized C-implementations (for both CPU and GPU) and will default to Python implementations. Performance will be severely degraded. To remove this warning, set Theano flags cxx to an empty string.\n", "WARNING (theano.configdefaults): install mkl with `conda install mkl-service`: No module named 'mkl'\n", "WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n" ] } ], "source": [ "import pandas as pd\n", "import gempy as gp\n", "import mplstereonet\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "id": "jewish-complex", "metadata": {}, "source": [ "\n", "# Data Preparation\n", "\n", "For this model, the only thing that needs to be done is loading the already created interface points and orientations. In the next tutorials, you will create the data yourself and process it further to make it usable for GemPy. \n", "\n", "\n", "## Importing Interface Points\n", "\n", "We are using the `pandas` library to load the interface points that were prepared beforehand and stored as CSV-file. The only information that is needed are the location of the interface point (`X`, `Y`, `Z`) and the `formation` it belongs to. You may have to adjust the `delimiter` when loading the file." ] }, { "cell_type": "code", "execution_count": 3, "id": "9c8c178c", "metadata": { "ExecuteTime": { "end_time": "2022-04-02T12:33:29.903687Z", "start_time": "2022-04-02T12:33:29.874543Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
XYZformation
0250250-250Layer1
1250500-250Layer1
2250750-250Layer1
3500250-50Layer1
4500500-50Layer1
\n", "
" ], "text/plain": [ " X Y Z formation\n", "0 250 250 -250 Layer1\n", "1 250 500 -250 Layer1\n", "2 250 750 -250 Layer1\n", "3 500 250 -50 Layer1\n", "4 500 500 -50 Layer1" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "interfaces = pd.read_csv('../../data/model2/model2_interfaces.csv', \n", " delimiter = ';')\n", "interfaces.head()" ] }, { "cell_type": "markdown", "id": "occupied-simulation", "metadata": {}, "source": [ "\n", "\n", "## Importing Orientations\n", "\n", "The orientations will also be loaded using `pandas`. In addition to the location and the formation the orientation belongs to, the dip value, azimuth value (dip direction) and a polarity value (mostly set to 1 by default) needs to be provided. As the model will feature folded layers, the dip and the azimuth are variable. Orientations are provided for all four modeled layers. " ] }, { "cell_type": "code", "execution_count": null, "id": "bd8aa29b", "metadata": { "ExecuteTime": { "end_time": "2022-04-02T12:33:29.935404Z", "start_time": "2022-04-02T12:33:29.906205Z" } }, "outputs": [], "source": [ "orientations = pd.read_csv('../../data/model2/model2_orientations.csv', \n", " delimiter=';')\n", "orientations.head()" ] }, { "cell_type": "markdown", "id": "70a49b50", "metadata": {}, "source": [ "\n", "\n", "### Visualizing the orientations using mplstereonet\n", "\n", "The open-source Python package `mplstereonet` ([mplstereonet](https://mplstereonet.readthedocs.io/en/latest/index.html)) is used to plot the orientations needed for constraining the model of the folded layers. The stereonet displays the planes and poles for the provided orientations as lower hemisphere plot. In this example, it is obvious that the layers dip with 45° towards the east and west, depending on the location of the orientation within the folded structure. \n", "\n", "There is also a [tutorial available for this task on the GemGIS Documentation page](https://gemgis.readthedocs.io/en/latest/getting_started/tutorial/17_plotting_orientations_with_mplstereonet.html)." ] }, { "cell_type": "code", "execution_count": null, "id": "d3dbad5c", "metadata": {}, "outputs": [], "source": [ "fig = plt.figure()\n", "ax = fig.add_subplot(111, projection='stereonet')\n", "strike, dip = orientations['azimuth'].values+90, orientations['dip'].values\n", "ax.plane(strike, dip,linewidth=2, color='black')\n", "ax.pole(strike, dip, color='black')\n", "ax.grid()\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "veterinary-acceptance", "metadata": {}, "source": [ "\n", "\n", "# GemPy Model Calculation\n", "\n", "The following part presents the main steps of creating a model in `GemPy`. \n", "\n", "The creation of a `GemPy` Model follows particular steps which will be performed in the following:\n", "\n", "1. Create new model: `gp.create_model()`\n", "2. Data Initiation: `gp.init_data()`\n", "3. Map Stack to Surfaces: `gp.map_stack_to_surfaces()`\n", "4. [...]\n", "5. Set the Interpolator: `gp.set_interpolator()`\n", "6. Computing the Model: `gp.compute_model()`\n", "\n", "\n", "\n", "## Creating the GemPy Model\n", "\n", "The first step is to create a new empty `GemPy` model by providing a name for it. " ] }, { "cell_type": "code", "execution_count": null, "id": "eb3ec438", "metadata": { "ExecuteTime": { "end_time": "2022-04-02T12:33:30.683572Z", "start_time": "2022-04-02T12:33:29.938400Z" } }, "outputs": [], "source": [ "geo_model = gp.create_model('Model2_Folded_Layers')\n", "geo_model" ] }, { "cell_type": "markdown", "id": "supreme-encoding", "metadata": {}, "source": [ "\n", "\n", "## Data Initiation\n", "\n", "During this step, the `extent` of the model (`xmin`, `xmax`, `ymin`, `ymax`, `zmin`, `zmax`) and the `resolution` in `X`, `Y`and `Z` direction (`res_x`, `res_y`, `res_z`, equal to the number of cells in each direction) will be set using lists of values. \n", "\n", "The interface points (`surface_points_df`) and orientations (`orientations_df`) will be passed as `pandas` `DataFrames`. " ] }, { "cell_type": "code", "execution_count": null, "id": "ae4bb825", "metadata": { "ExecuteTime": { "end_time": "2022-04-02T12:33:30.936748Z", "start_time": "2022-04-02T12:33:30.685569Z" } }, "outputs": [], "source": [ "gp.init_data(geo_model=geo_model, \n", " extent=[0, 1000, 0, 1000, -600, 0], \n", " resolution=[100, 100, 100],\n", " surface_points_df=interfaces,\n", " orientations_df=orientations,\n", " default_values=True)" ] }, { "cell_type": "markdown", "id": "armed-story", "metadata": {}, "source": [ "\n", "\n", "## Inspecting the Surfaces\n", "\n", "The model consists of four different layers or surfaces now which all belong to the `Default series`. During the next step, the proper `Series` will be assigned to the surfaces. Using the `surfaces`-attribute again, we can check which layers were loaded." ] }, { "cell_type": "code", "execution_count": null, "id": "9e8519ea", "metadata": { "ExecuteTime": { "end_time": "2022-04-02T12:33:31.030680Z", "start_time": "2022-04-02T12:33:30.939798Z" } }, "outputs": [], "source": [ "geo_model.surfaces" ] }, { "cell_type": "markdown", "id": "980fdbc6", "metadata": {}, "source": [ "\n", "\n", "## Inspecting the Input Data\n", "\n", "The loaded interface points and orientations can again be inspected using the `surface_points`- and `orientations`-attributes. Using the `df`-attribute of this object will convert the displayed table in a `pandas` `DataFrame`." ] }, { "cell_type": "code", "execution_count": null, "id": "07fc390c", "metadata": {}, "outputs": [], "source": [ "geo_model.surface_points.df.head()" ] }, { "cell_type": "code", "execution_count": null, "id": "9a171ee0", "metadata": {}, "outputs": [], "source": [ "geo_model.orientations.df.head()" ] }, { "cell_type": "markdown", "id": "occupational-coaching", "metadata": {}, "source": [ "\n", "\n", "## Map Stack to Surfaces\n", "\n", "During this step, all four layers of the model are assigned to the `Strata1` series. We know that the layers modeled here are parallel. If the layers were not parallel as shown in the next models, multiple series would be defined. We will also add a `Basement` here (`geo_model.add_surfaces('Basement')`). The order within one series also defines the age relations within this series and has to be according to the depositional events of the layers." ] }, { "cell_type": "code", "execution_count": null, "id": "ca6dfa3e", "metadata": { "ExecuteTime": { "end_time": "2022-04-02T12:33:31.174730Z", "start_time": "2022-04-02T12:33:31.033684Z" } }, "outputs": [], "source": [ "gp.map_stack_to_surfaces(geo_model,\n", " {\n", " 'Strata1': ('Layer1', 'Layer2', 'Layer3', 'Layer4'), \n", " },\n", " remove_unused_series=True)\n", "geo_model.add_surfaces('Basement')\n", "geo_model.surfaces" ] }, { "cell_type": "code", "execution_count": null, "id": "3f6e7b7f", "metadata": {}, "outputs": [], "source": [ "geo_model.stack" ] }, { "cell_type": "markdown", "id": "directed-immigration", "metadata": {}, "source": [ "\n", "\n", "## Plotting the input data in 2D using Matplotlib\n", "\n", "The input data can now be visualized in 2D using `matplotlib`. This might for example be useful to check if all points and measurements are defined the way we want them to. Using the function `plot_2d()`, we attain a 2D projection of our data points onto a plane of chosen direction (we can choose this attribute to be either `'x'`, `'y'`, or `'z'`)." ] }, { "cell_type": "code", "execution_count": null, "id": "ee36f42b", "metadata": { "ExecuteTime": { "end_time": "2022-04-02T12:33:31.493549Z", "start_time": "2022-04-02T12:33:31.178718Z" } }, "outputs": [], "source": [ "gp.plot_2d(geo_model, \n", " direction='z', \n", " show_lith=False, \n", " show_boundaries=False)\n", "plt.grid()" ] }, { "cell_type": "markdown", "id": "objective-standard", "metadata": {}, "source": [ "\n", "\n", "## Plotting the input data in 3D using PyVista\n", "\n", "The input data can also be viszualized using the `pyvista` package. In this view, the interface points are visible as well as the orientations (marked as arrows) which indicate the normals of each orientation value. " ] }, { "cell_type": "code", "execution_count": null, "id": "03cfc42d", "metadata": { "ExecuteTime": { "end_time": "2022-04-02T12:33:32.096033Z", "start_time": "2022-04-02T12:33:31.497542Z" } }, "outputs": [], "source": [ "gp.plot_3d(geo_model, \n", " image=False, \n", " plotter_type='basic', \n", " notebook=True)" ] }, { "cell_type": "markdown", "id": "elementary-camping", "metadata": {}, "source": [ "\n", "## Setting the interpolator\n", "\n", "Once we have made sure that we have defined all our primary information, we can continue with the next step towards creating our geological model: preparing the input data for interpolation.\n", "\n", "Setting the interpolator is necessary before computing the actual model. Here, the most important kriging parameters can be defined. " ] }, { "cell_type": "code", "execution_count": null, "id": "8e27b70d", "metadata": { "ExecuteTime": { "end_time": "2022-04-02T12:33:35.201651Z", "start_time": "2022-04-02T12:33:32.098033Z" } }, "outputs": [], "source": [ "gp.set_interpolator(geo_model,\n", " compile_theano=True,\n", " theano_optimizer='fast_compile',\n", " verbose=[],\n", " update_kriging=False\n", " )" ] }, { "cell_type": "markdown", "id": "innovative-cradle", "metadata": {}, "source": [ "\n", "\n", "## Computing the model\n", "\n", "At this point, we have all we need to compute our full model via `gp.compute_model()`. By default, this will return two separate solutions in the form of arrays. The first provides information on the lithological formations, the second on the fault network in the model, which is not present in this example. " ] }, { "cell_type": "code", "execution_count": null, "id": "5b7b2377", "metadata": { "ExecuteTime": { "end_time": "2022-04-02T12:34:25.093460Z", "start_time": "2022-04-02T12:33:35.205651Z" } }, "outputs": [], "source": [ "sol = gp.compute_model(geo_model, \n", " compute_mesh=True)" ] }, { "cell_type": "code", "execution_count": null, "id": "d161bace", "metadata": {}, "outputs": [], "source": [ "sol" ] }, { "cell_type": "code", "execution_count": null, "id": "a03ac213", "metadata": {}, "outputs": [], "source": [ "geo_model.solutions" ] }, { "cell_type": "markdown", "id": "protected-parcel", "metadata": {}, "source": [ "\n", "\n", "# Model Visualization and Post-Processing\n", "\n", "\n", "\n", "## Visulazing Cross Sections of the computed model\n", "\n", "Cross sections in different `direction`s and at different `cell_number`s can be displayed. Here, we see the folded layers of the model in the different directions. " ] }, { "cell_type": "code", "execution_count": null, "id": "eb54b87a", "metadata": { "ExecuteTime": { "end_time": "2022-04-02T12:34:26.391340Z", "start_time": "2022-04-02T12:34:25.096457Z" } }, "outputs": [], "source": [ "gp.plot_2d(geo_model, \n", " direction=['x', 'x', 'y', 'y'], \n", " cell_number=[25, 50, 25, 75], \n", " show_topography=False, \n", " show_data=True)" ] }, { "cell_type": "markdown", "id": "bd7ca4a2", "metadata": {}, "source": [ "Next to the lithology data, we can also plot the calculated scalar field." ] }, { "cell_type": "code", "execution_count": null, "id": "ad70273d", "metadata": {}, "outputs": [], "source": [ "gp.plot_2d(geo_model, direction='y', show_data=False, show_scalar=True, show_lith=False)" ] }, { "cell_type": "markdown", "id": "several-artwork", "metadata": {}, "source": [ "\n", "\n", "# Visualizing the computed model in 3D\n", "\n", "The computed model can be visualized in 3D using the `pyvista` library. Setting `notebook=False` will open an interactive windows and the model can be rotated and zooming is possible. " ] }, { "cell_type": "code", "execution_count": null, "id": "838f0e34", "metadata": { "ExecuteTime": { "end_time": "2022-04-02T12:34:27.136575Z", "start_time": "2022-04-02T12:34:26.394341Z" } }, "outputs": [], "source": [ "gpv = gp.plot_3d(geo_model, \n", " image=False, \n", " show_topography=True,\n", " plotter_type='basic', \n", " notebook=True, \n", " show_lith=True)" ] }, { "cell_type": "markdown", "id": "ee08db71", "metadata": {}, "source": [ "\n", "# Conclusions\n", "\n", "
\n", "In this tutorial, you have learnt the following:
\n", "- Get a better undestanding of what orientations mean in GemPy and how to plot them using mplstereonet
\n", "- How to build a simple model consisting of folded layers belonging to one Series
\n", "\n", "
\n", "\n", "\n", "\n", "# Outlook\n", "\n", "
\n", "In the next tutorial, you will learn the following:
\n", "- Get an understanding of how faults are modeled/displayed in GemPy
\n", "- How to build a simple model consisting of faulted layers belonging to one Series
\n", "\n", "\n", "
\n", "\n", "[Take me to the next notebook on Github](https://nbviewer.org/github/cgre-aachen/gemgis_data/blob/main/notebooks/01_basic_modeling/model3_Faulted_Layers.ipynb)\n", "\n", "[Take me to the next notebook locally](model3_Faulted_Layers.ipynb)\n", "\n", "\n", "\n", "\n", "\n", "## Licensing\n", "\n", "Institute for Computational Geoscience, Geothermics and Reservoir Geophysics, RWTH Aachen University & Fraunhofer IEG, Fraunhofer Research Institution for Energy Infrastructures and Geothermal Systems IEG, Authors: Alexander Juestel. For more information contact: alexander.juestel(at)ieg.fraunhofer.de\n", "\n", "All notebooks are licensed under a Creative Commons Attribution 4.0 International License (CC BY 4.0, http://creativecommons.org/licenses/by/4.0/). References for each displayed map are provided. Most of the maps originate from the books of [Powell (1992)](https://link.springer.com/book/9783540586074) and [Bennison (1990)](https://link.springer.com/book/10.1007/978-1-4615-9630-1). References for maps with unknown origin will gladly be added." ] } ], "metadata": { "hide_input": false, "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.15" }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 5 }