{ "cells": [ { "cell_type": "markdown", "id": "arranged-tenant", "metadata": {}, "source": [ "\n", "\n", "#
From Maps to Models - Tutorials for structural geological modeling using GemPy and GemGIS
\n", "\n", "# Example 2 - Folded Layers\n", "\n", "This example will show how to convert the geological map below using ``GemGIS`` to a `GemPy` model. This example is based on digitized data. The area is 4016 m wide (W-E extent) and 2790 m high (N-S extent). The vertical model extent varies from -250 m to 1500 m. The model represents folded layers (blue to light green) above a basement unit (light blue).\n", "\n", "
\n", "In this tutorial, you will learn the following:
\n", "- How to build your fifth GemPy model with input data generated through GemGIS
\n", "- How to sample rasters (interfaces/orientations)
\n", "- How to calculate raster properties and how to perform raster operations
\n", "- How to visualize the input data in PyVista
\n", "
\n", "\n", "## Your Tasks\n", "\n", "\n", "1. Georeference the map in QGIS given the dimensions above using the coordinate reference system with the EPSG code 4326\n", "2. Digitize the layer boundaries (including a `formation` column) and the topographic lines (including a `Z` column)\n", "3. Digitize so-called strike lines for the different layers. The orientations used for `GemPy` will be calculated from the strike lines.\n", "\n", "## Contents\n", "\n", "1. [Installing GemPy and GemGIS](#installing-gempy)\n", "2. [Importing Libraries](#importing-libraries)\n", "3. [Data Preparation](#data-preparation)\n", " 1. [Creating Digital Elevation Model from Contour lines](#creating-digital-elevation-model-from-contour-lines)\n", " 1. [Loading Contour Lines](#loading-contour-lines)\n", " 2. [Plotting Contour Lines](#plotting-contour-lines)\n", " 3. [Interpolating Contour Lines](#interpolating-contour-lines)\n", " 4. [Plotting the raster](#plotting-the-raster)\n", " 5. [Saving the raster to disc](#saving-the-raster-to-disc)\n", " 6. [Opening raster](#opening-raster)\n", " 2. [Processing Stratigraphic Boundaries](#processing-stratigraphic-boundaries)\n", " 1. [Opening Stratigraphic Boundaries](#opening-stratigraphic-boundaries)\n", " 2. [Plotting Stratigraphic Boundaries](#plotting-stratigraphic-boundaries)\n", " 3. [Extracting Z coordinates from Digital Elevation Model](#extracting-z-coordinates-from-digital-elevation-model)\n", " 4. [Plotting the Interface Points](#plotting-the-interface-points)\n", " 3. [Processing Orientations](#processing-orientations)\n", " 1. [Orientations from Strike Lines](#orientations-from-strike-lines)\n", " 2. [Loading Strike Lines](#loading-strike-lines) \n", " 2. [Plotting Strike Lines](#plotting-strike-lines) \n", " 4. [Calculate Orientations for each formation](#calculate-orientations-for-each-formation)\n", " 5. [Merging Orientations for GemPy](#merging-orientations-for-gempy)\n", " 6. [Plotting Orientations](#plotting-the-orientations)\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. [Showing the Number of Data Points](#showing-the-number-of-data-points)\n", " 7. [Loading Digital Elevation Model](#loading-digital-elevation-model)\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 the computed model in 2D](#visualizing-the-computed-model-in-2d)\n", " 2. [Visualizing the computed model in 3D](#visualizing-the-computed-model-in-3d) \n", " 3. [Post Processing](#post-processing)\n", " 1. [Sampling Interfaces and Orientations from Raster](#sampling-raster)\n", " 2. [Calculating Raster Properties](#calculating-raster-properties)\n", " 3. [Performing Raster Operations](#performing-raster-operations)\n", " 4. [Visualizing Input Data with PyVista](#visualizing-input-data-with-pyvista)\n", " \n", "\n", "Source: Bennison, G.M. (1988): An Introduction to Geological Structures and Maps, page 26, figure 11, Springer Verlag Berlin, Heidelberg, New York, ISBN: 978-1-4615-9632-5" ] }, { "cell_type": "markdown", "id": "alike-configuration", "metadata": {}, "source": [ "\n", "\n", "# Installing GemPy and GemGIS\n", "\n", "If you have not installed `GemPy` yet, please follow the [GemPy installation instructions](https://docs.gempy.org/installation.html) and the [GemGIS installation instructions](https://gemgis.readthedocs.io/en/latest/getting_started/installation.html). If you encounter any issues, feel free to open a new discussion at [GemPy Discussions](https://github.com/cgre-aachen/gempy/discussions) or [GemGIS Discussions](https://github.com/cgre-aachen/gemgis/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) or [GemGIS Issues](https://github.com/cgre-aachen/gemgis/issues). There, the `GemPy` and `GemGIS` development teams will help you out. " ] }, { "cell_type": "markdown", "id": "allied-vessel", "metadata": {}, "source": [ "\n", "\n", "# Importing Libraries\n", "\n", "For this notebook, we need the `geopandas` library for the data preparation, `rasterio` for dealing with the created digital elevation model, `matplotlib` for plotting, `numpy` for some numerical calculations, `pandas` for manipulating `DataFrames` and of course the `gempy` and `gemgis` libraries. Any warnings that may appear can be ignored for now. The file path is set to load the data provided for this tutorial. " ] }, { "cell_type": "code", "execution_count": null, "id": "polar-appliance", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:23.414255Z", "start_time": "2022-04-06T09:52:23.401235Z" } }, "outputs": [], "source": [ "import geopandas as gpd\n", "import rasterio\n", "import warnings\n", "warnings.filterwarnings(\"ignore\")\n", "import gemgis as gg\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import gempy as gp\n", "import pyvista as pv\n", "import pandas as pd" ] }, { "cell_type": "code", "execution_count": null, "id": "healthy-queens", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:23.430256Z", "start_time": "2022-04-06T09:52:23.416237Z" } }, "outputs": [], "source": [ "file_path = '../../data/example07_folded_layers/'" ] }, { "cell_type": "markdown", "id": "impressed-database", "metadata": {}, "source": [ "\n", "# Data Preparation\n", "\n", "At his point, you should have the topographic contour lines (including a `Z` column) and the layer boundaries (including a `formation` column) digitized. If not, please generate the data before continuing with this tutorial. \n", "\n", "\n", "## Creating Digital Elevation Model from Contour Lines\n", "\n", "The digital elevation model (DEM) will be created by interpolating the contour lines digitized from the georeferenced map using the `SciPy` [Radial Basis Function interpolation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.Rbf.html) wrapped in `GemGIS`. The respective function used for that is `gg.vector.interpolate_raster()`. \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/05_interpolating_rasters.html).\n", "\n", "\n", "\n", "Source: Bennison, G.M. (1988): An Introduction to Geological Structures and Maps, page 26, figure 11, Springer Verlag Berlin, Heidelberg, New York, ISBN: 978-1-4615-9632-5\n", "\n", "\n", "\n", "### Loading contour lines\n", "\n", "First, the contour lines are loaded using `GeoPandas`. Please provide here the name of your shape file containing the digitized topographic contour lines. " ] }, { "cell_type": "code", "execution_count": null, "id": "informed-tomorrow", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:23.699300Z", "start_time": "2022-04-06T09:52:23.622282Z" } }, "outputs": [], "source": [ "topo = gpd.read_file(file_path + 'topo7.shp')\n", "topo.head()" ] }, { "cell_type": "markdown", "id": "14e418ef", "metadata": {}, "source": [ "\n", "\n", "### Plotting the contour lines\n", "\n", "The contour lines are plotted using the built-in plotting function of `GeoPandas`. " ] }, { "cell_type": "code", "execution_count": null, "id": "bcc57724", "metadata": {}, "outputs": [], "source": [ "topo.plot(column='Z', aspect=1, legend=True, cmap='gist_earth')" ] }, { "cell_type": "markdown", "id": "inclusive-florence", "metadata": {}, "source": [ "\n", "\n", "### Interpolating the contour lines\n", "\n", "The digital elevation model (DEM) will be created by interpolating the contour lines digitized from the georeferenced map using the `SciPy` [Radial Basis Function interpolation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.Rbf.html) wrapped in `GemGIS`. The respective function used for that is `gg.vector.interpolate_raster()`. \n" ] }, { "cell_type": "code", "execution_count": null, "id": "revolutionary-calibration", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:24.865684Z", "start_time": "2022-04-06T09:52:23.702304Z" } }, "outputs": [], "source": [ "topo_raster = gg.vector.interpolate_raster(gdf=topo, value='Z', method='rbf', res=10)" ] }, { "cell_type": "markdown", "id": "experimental-press", "metadata": {}, "source": [ "\n", "\n", "### Plotting the raster\n", "\n", "The interpolated digital elevation model can be displayed using `matplotlib` and its `plt.imshow()` function and by providing the extent of the raster to align it with the contour lines. " ] }, { "cell_type": "code", "execution_count": null, "id": "frozen-gazette", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:25.184063Z", "start_time": "2022-04-06T09:52:24.867686Z" } }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "\n", "from mpl_toolkits.axes_grid1 import make_axes_locatable\n", "fix, ax = plt.subplots(1, figsize=(5,5))\n", "topo.plot(ax=ax, aspect='equal', column='Z', cmap='gist_earth')\n", "im = plt.imshow(topo_raster, origin='lower', extent=[0, 4016, 0, 2790], cmap='gist_earth')\n", "divider = make_axes_locatable(ax)\n", "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n", "cbar = plt.colorbar(im, cax=cax)\n", "cbar.set_label('Altitude [m]')\n", "ax.set_xlabel('X [m]')\n", "ax.set_ylabel('Y [m]')\n", "ax.set_xlim(0, 4016)\n", "ax.set_ylim(0, 2790)" ] }, { "cell_type": "markdown", "id": "common-franklin", "metadata": {}, "source": [ "\n", "\n", "### Saving the raster to disc\n", "\n", "After the interpolation of the contour lines, the raster is saved to disc using `gg.raster.save_as_tiff()`. The function will not be executed as a raster is already provided with the example data. " ] }, { "cell_type": "raw", "id": "renewable-walker", "metadata": { "ExecuteTime": { "end_time": "2021-03-29T15:10:34.847704Z", "start_time": "2021-03-29T15:10:34.720676Z" } }, "source": [ "gg.raster.save_as_tiff(raster=topo_raster, path=file_path + 'raster7.tif', extent=[0, 4016, 0, 2790], crs='EPSG:4326', overwrite_file=True)" ] }, { "cell_type": "markdown", "id": "under-pursuit", "metadata": {}, "source": [ "\n", "\n", "### Opening Raster\n", "\n", "The previously computed and saved raster can now be opened using rasterio. " ] }, { "cell_type": "code", "execution_count": null, "id": "silver-defendant", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:25.232192Z", "start_time": "2022-04-06T09:52:25.186063Z" } }, "outputs": [], "source": [ "topo_raster = rasterio.open(file_path + 'raster7.tif')" ] }, { "cell_type": "markdown", "id": "central-option", "metadata": { "ExecuteTime": { "end_time": "2021-03-28T12:43:27.113805Z", "start_time": "2021-03-28T12:43:27.067793Z" } }, "source": [ "\n", "\n", "## Processing Stratigraphic Boundaries\n", "\n", "The interface points will be extracted from LineStrings digitized from the georeferenced map using QGIS. It is important to provide a `formation` name for each layer boundary. Up until now, only the `X` and `Y` position are stored in the vertices of the LineStrings. Using the digital elevation model created already, we will now sample the elevation model at the locations of the vertices to extract the height at this point as the stratigraphic boundary was mapped at the surface.\n", "\n", "\n", "Source: Bennison, G.M. (1988): An Introduction to Geological Structures and Maps, page 26, figure 11, Springer Verlag Berlin, Heidelberg, New York, ISBN: 978-1-4615-9632-5" ] }, { "cell_type": "markdown", "id": "6cb95b25", "metadata": {}, "source": [ "\n", "\n", "### Opening Stratigraphic Boundaries\n", "\n", "The stratigraphic units are opened using `GeoPandas`." ] }, { "cell_type": "code", "execution_count": null, "id": "activated-namibia", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:25.472248Z", "start_time": "2022-04-06T09:52:25.426221Z" } }, "outputs": [], "source": [ "interfaces = gpd.read_file(file_path + 'interfaces7.shp')\n", "interfaces.head()" ] }, { "cell_type": "markdown", "id": "545cc9d8", "metadata": {}, "source": [ "\n", "\n", "### Plotting Stratigraphic Boundaries" ] }, { "cell_type": "code", "execution_count": null, "id": "c7fdd21d", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(1, figsize=(5,5))\n", "\n", "interfaces.plot(ax=ax, column='formation', legend=True, aspect='equal')\n", "\n", "plt.grid()\n", "ax.set_xlabel('X [m]')\n", "ax.set_ylabel('Y [m]')" ] }, { "cell_type": "markdown", "id": "prospective-drink", "metadata": {}, "source": [ "\n", "\n", "### Extracting Z coordinates from Digital Elevation Model\n", "\n", "The vertical position of the interface point will be extracted from the digital elevation model using the `GemGIS` function `gg.vector.extract_xyz()`. The resulting GeoDataFrame now contains single points including the information about the respective `formation` as well as the `X`, `Y`, and `Z` location. This is all we need as preparational steps to generate input data for `GemPy`. \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/02_extract_xyz.html)." ] }, { "cell_type": "code", "execution_count": null, "id": "stopped-saint", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:25.583271Z", "start_time": "2022-04-06T09:52:25.474248Z" } }, "outputs": [], "source": [ "interfaces_coords = gg.vector.extract_xyz(gdf=interfaces, dem=topo_raster)\n", "interfaces_coords = interfaces_coords.sort_values(by='formation', ascending=False)\n", "interfaces_coords.head()" ] }, { "cell_type": "markdown", "id": "valued-coast", "metadata": {}, "source": [ "\n", "\n", "### Plotting the Interface Points\n", "\n", "The interface points incuding their altitude (Z-) values and the digitized LineString can be plotted using `matplotlib`. " ] }, { "cell_type": "code", "execution_count": null, "id": "involved-association", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:25.950339Z", "start_time": "2022-04-06T09:52:25.585264Z" } }, "outputs": [], "source": [ "fig, ax = plt.subplots(1, figsize=(5,5))\n", "\n", "interfaces.plot(ax=ax, column='formation', legend=True, aspect='equal')\n", "interfaces_coords.plot(ax=ax, column='formation', legend=True, aspect='equal')\n", "plt.grid()\n", "ax.set_xlabel('X [m]')\n", "ax.set_ylabel('Y [m]')\n", "ax.set_xlim(0, 4016)\n", "ax.set_ylim(0, 2790)" ] }, { "cell_type": "markdown", "id": "inner-mathematics", "metadata": {}, "source": [ "\n", "\n", "## Processing Orientations\n", "\n", "For this example, orientations must be calculated yourself. They will be calculated using functions implemented in GemGIS and the previously digitized strike lines. \n", "\n", "\n", "Source: Bennison, G.M. (1988): An Introduction to Geological Structures and Maps, page 26, figure 11, Springer Verlag Berlin, Heidelberg, New York, ISBN: 978-1-4615-9632-5\n", "\n", "\n", "### Orientations from Strike Lines\n", "\n", "Strike lines connect outcropping stratigraphic boundaries (interfaces) of the same altitude. In other words: the intersections between topographic contours and stratigraphic boundaries at the surface. The height difference and the horizontal difference between two digitized lines is used to calculate the dip and azimuth and hence an orientation that is necessary for `GemPy`.\n", "\n", "The calculation of orientations from strike lines has been implemented into `GemPy` for simple cases like these. In order to calculate the orientations, each set of strikes lines/LineStrings for one formation must be given an id number next to the altitude of the strike line. The id field is already predefined in QGIS. The strike line with the lowest altitude gets the id number `1`, the strike line with the highest altitude the the number according to the number of digitized strike lines. It is currently recommended to use one set of strike lines for each structural element of one formation as illustrated. \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", "\n", "Source: Powell, D. (1995): Interpretation geologischer Strukturen durch Karten - Eine praktische Anleitung mit Aufgaben und Lösungen, page 14, figure 8, Springer Verlag Berlin, Heidelberg, New York, ISBN: 978-3-540-58607-4." ] }, { "cell_type": "markdown", "id": "8b1e7bd3", "metadata": {}, "source": [ "\n", "### Loading Strike Lines\n", "\n", "We are using `GeoPandas` to load the strike lines. " ] }, { "cell_type": "code", "execution_count": null, "id": "configured-evaluation", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:26.189953Z", "start_time": "2022-04-06T09:52:26.143939Z" }, "scrolled": true }, "outputs": [], "source": [ "strikes = gpd.read_file(file_path + 'strikes7.shp')\n", "strikes.head()" ] }, { "cell_type": "markdown", "id": "ecec2a31", "metadata": {}, "source": [ "\n", "\n", "### Plotting Strike Lines\n", "\n", "The strike lines can be plotted using `matplotlib`. " ] }, { "cell_type": "code", "execution_count": null, "id": "633b3bbb", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(1, figsize=(5,5))\n", "\n", "strikes.plot(ax=ax,column='id', aspect=1)\n", "interfaces.plot(ax=ax, column='formation', legend=True, aspect='equal')\n", "ax.set_xlabel('X [m]')\n", "ax.set_ylabel('Y [m]')\n", "ax.grid()" ] }, { "cell_type": "markdown", "id": "imposed-elite", "metadata": {}, "source": [ "\n", "\n", "### Calculate Orientations for each formation\n", "\n", "The calculations will be calculated using the function `gg.vector.calculate_orientations_from_strike_lines()` where the strike lines for each single formation will be provided and calculated separately. The result is a `GeoDataFrame` ready to be used in `GemPy`." ] }, { "cell_type": "code", "execution_count": null, "id": "durable-fault", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:26.269983Z", "start_time": "2022-04-06T09:52:26.191951Z" } }, "outputs": [], "source": [ "orientations_b = gg.vector.calculate_orientations_from_strike_lines(gdf=strikes[strikes['formation'] == 'B'].sort_values(by='Z', ascending=True).reset_index())\n", "orientations_b" ] }, { "cell_type": "code", "execution_count": null, "id": "tropical-ticket", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:26.350002Z", "start_time": "2022-04-06T09:52:26.271971Z" } }, "outputs": [], "source": [ "orientations_b1 = gg.vector.calculate_orientations_from_strike_lines(gdf=strikes[strikes['formation'] == 'B1'].sort_values(by='Z', ascending=True).reset_index())\n", "orientations_b1" ] }, { "cell_type": "code", "execution_count": null, "id": "sweet-nation", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:26.414013Z", "start_time": "2022-04-06T09:52:26.352004Z" } }, "outputs": [], "source": [ "orientations_c = gg.vector.calculate_orientations_from_strike_lines(gdf=strikes[strikes['formation'] == 'C'].sort_values(by='Z', ascending=True).reset_index())\n", "orientations_c" ] }, { "cell_type": "code", "execution_count": null, "id": "bridal-nitrogen", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:26.493030Z", "start_time": "2022-04-06T09:52:26.416003Z" } }, "outputs": [], "source": [ "orientations_c1 = gg.vector.calculate_orientations_from_strike_lines(gdf=strikes[strikes['formation'] == 'C1'].sort_values(by='Z', ascending=True).reset_index())\n", "orientations_c1" ] }, { "cell_type": "code", "execution_count": null, "id": "coordinate-topic", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:26.573056Z", "start_time": "2022-04-06T09:52:26.495030Z" } }, "outputs": [], "source": [ "orientations_d = gg.vector.calculate_orientations_from_strike_lines(gdf=strikes[strikes['formation'] == 'D'].sort_values(by='Z', ascending=True).reset_index())\n", "orientations_d" ] }, { "cell_type": "code", "execution_count": null, "id": "healthy-refund", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:26.637071Z", "start_time": "2022-04-06T09:52:26.575046Z" } }, "outputs": [], "source": [ "orientations_d1 = gg.vector.calculate_orientations_from_strike_lines(gdf=strikes[strikes['formation'] == 'D1'].sort_values(by='Z', ascending=True).reset_index())\n", "orientations_d1" ] }, { "cell_type": "code", "execution_count": null, "id": "neural-operations", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:26.716089Z", "start_time": "2022-04-06T09:52:26.642055Z" } }, "outputs": [], "source": [ "orientations_e = gg.vector.calculate_orientations_from_strike_lines(gdf=strikes[strikes['formation'] == 'E'].sort_values(by='Z', ascending=True).reset_index())\n", "orientations_e" ] }, { "cell_type": "code", "execution_count": null, "id": "pleasant-appearance", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:26.796107Z", "start_time": "2022-04-06T09:52:26.720070Z" } }, "outputs": [], "source": [ "orientations_f = gg.vector.calculate_orientations_from_strike_lines(gdf=strikes[strikes['formation'] == 'F'].sort_values(by='Z', ascending=True).reset_index())\n", "orientations_f" ] }, { "cell_type": "code", "execution_count": null, "id": "medical-hartford", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:26.860111Z", "start_time": "2022-04-06T09:52:26.798093Z" } }, "outputs": [], "source": [ "orientations_f1 = gg.vector.calculate_orientations_from_strike_lines(gdf=strikes[strikes['formation'] == 'F1'].sort_values(by='Z', ascending=True).reset_index())\n", "orientations_f1" ] }, { "cell_type": "code", "execution_count": null, "id": "developing-probability", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:26.939140Z", "start_time": "2022-04-06T09:52:26.862107Z" } }, "outputs": [], "source": [ "orientations_g = gg.vector.calculate_orientations_from_strike_lines(gdf=strikes[strikes['formation'] == 'G'].sort_values(by='Z', ascending=True).reset_index())\n", "orientations_g" ] }, { "cell_type": "code", "execution_count": null, "id": "white-bumper", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:27.018159Z", "start_time": "2022-04-06T09:52:26.941140Z" } }, "outputs": [], "source": [ "orientations_g1 = gg.vector.calculate_orientations_from_strike_lines(gdf=strikes[strikes['formation'] == 'G1'].sort_values(by='Z', ascending=True).reset_index())\n", "orientations_g1" ] }, { "cell_type": "code", "execution_count": null, "id": "faced-spelling", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:27.098177Z", "start_time": "2022-04-06T09:52:27.020144Z" } }, "outputs": [], "source": [ "orientations_h = gg.vector.calculate_orientations_from_strike_lines(gdf=strikes[strikes['formation'] == 'H'].sort_values(by='Z', ascending=True).reset_index())\n", "orientations_h" ] }, { "cell_type": "code", "execution_count": null, "id": "fitted-produce", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:27.162191Z", "start_time": "2022-04-06T09:52:27.100176Z" } }, "outputs": [], "source": [ "orientations_h1 = gg.vector.calculate_orientations_from_strike_lines(gdf=strikes[strikes['formation'] == 'H1'].sort_values(by='Z', ascending=True).reset_index())\n", "orientations_h1" ] }, { "cell_type": "markdown", "id": "friendly-token", "metadata": {}, "source": [ "\n", "\n", "### Merging Orientations for GemPy\n", "\n", "Since `GemPy` only takes one `DataFrame` for the necessary orientations, the single `DataFrames` are concatenated using `pd.concat()`. " ] }, { "cell_type": "code", "execution_count": null, "id": "sitting-reproduction", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:27.194187Z", "start_time": "2022-04-06T09:52:27.164174Z" } }, "outputs": [], "source": [ "import pandas as pd\n", "orientations = pd.concat([orientations_b, orientations_b1, orientations_c, orientations_c1, orientations_d, orientations_d1, orientations_e, orientations_f, orientations_f1, orientations_g, orientations_g1, orientations_h, orientations_h1]).reset_index()\n", "orientations['formation'] = ['B', 'B', 'B', 'B', 'B', 'B', 'B', 'C', 'C', 'C', 'C', 'C', 'C', 'D', 'D', 'D', 'D', 'D',\n", " 'E', 'E', 'E', 'F', 'F', 'F', 'F', 'F', 'G', 'G', 'G', 'G', 'G', 'G', 'G', 'H', 'H', 'H', 'H']\n", "orientations.head()" ] }, { "cell_type": "markdown", "id": "textile-christianity", "metadata": {}, "source": [ "\n", "\n", "### Plotting the Orientations\n", "\n", "The orientations can be plotted using `matplotlib`." ] }, { "cell_type": "code", "execution_count": null, "id": "detected-certificate", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:27.652768Z", "start_time": "2022-04-06T09:52:27.196183Z" } }, "outputs": [], "source": [ "fig, ax = plt.subplots(1, figsize=(5,5))\n", "\n", "interfaces.plot(ax=ax, column='formation', legend=True, aspect='equal')\n", "interfaces_coords.plot(ax=ax, column='formation', legend=True, aspect='equal')\n", "orientations.plot(ax=ax, color='red', aspect='equal')\n", "plt.grid()\n", "ax.set_xlabel('X [m]')\n", "ax.set_ylabel('Y [m]')\n", "ax.set_xlim(0, 4016)\n", "ax.set_ylim(0, 2790)" ] }, { "cell_type": "markdown", "id": "wired-vietnamese", "metadata": {}, "source": [ "\n", "\n", "# GemPy Model Calculation\n", "\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": "binding-fishing", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:30.279842Z", "start_time": "2022-04-06T09:52:29.387461Z" } }, "outputs": [], "source": [ "geo_model = gp.create_model('Model7')\n", "geo_model" ] }, { "cell_type": "markdown", "id": "stopped-poster", "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": "assisted-nothing", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:30.501885Z", "start_time": "2022-04-06T09:52:30.281843Z" } }, "outputs": [], "source": [ "gp.init_data(geo_model, [0, 4016, 0, 2790, -250, 1500], [100, 100, 100],\n", " surface_points_df=interfaces_coords[interfaces_coords['Z']!=0],\n", " orientations_df=orientations,\n", " default_values=True)" ] }, { "cell_type": "markdown", "id": "integrated-pension", "metadata": {}, "source": [ "\n", "\n", "## Inspecting the Surfaces\n", "\n", "The model consists of seven 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": "cutting-energy", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:30.596984Z", "start_time": "2022-04-06T09:52:30.502894Z" } }, "outputs": [], "source": [ "geo_model.surfaces" ] }, { "cell_type": "markdown", "id": "9fffbc9a", "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": "5eb18c24", "metadata": {}, "outputs": [], "source": [ "geo_model.surface_points.df.head()" ] }, { "cell_type": "code", "execution_count": null, "id": "d1fb46c3", "metadata": {}, "outputs": [], "source": [ "geo_model.orientations.df.head()" ] }, { "cell_type": "markdown", "id": "younger-young", "metadata": {}, "source": [ "\n", "\n", "## Map Stack to Surfaces\n", "\n", "During this step, all seven 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": "formed-thought", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:30.709027Z", "start_time": "2022-04-06T09:52:30.598994Z" } }, "outputs": [], "source": [ "gp.map_stack_to_surfaces(geo_model,\n", " { \n", " 'Strata1': ('H', 'G', 'F', 'E', 'D', 'C', 'B')\n", " },\n", " remove_unused_series=True)\n", "geo_model.add_surfaces('A')" ] }, { "cell_type": "markdown", "id": "subtle-techno", "metadata": {}, "source": [ "\n", "\n", "## Showing the Number of Data Points\n", "\n", "You can also return the number of interfaces and orientations for each formation using `gg.utils.show_number_of_data_points()`" ] }, { "cell_type": "code", "execution_count": null, "id": "forty-treaty", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:30.772042Z", "start_time": "2022-04-06T09:52:30.711010Z" } }, "outputs": [], "source": [ "gg.utils.show_number_of_data_points(geo_model=geo_model)" ] }, { "cell_type": "markdown", "id": "medieval-deployment", "metadata": {}, "source": [ "\n", "\n", "## Loading Digital Elevation Model\n", "\n", "`GemPy` is capable of including a topography into the modeling process. Here, we use the topography that we have interpolated in one of the previous steps. `GemPy` takes the file path of the raster/digital elevation model and loads it as grid into the `geo_model` object." ] }, { "cell_type": "code", "execution_count": null, "id": "subsequent-technique", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:31.185919Z", "start_time": "2022-04-06T09:52:30.774021Z" } }, "outputs": [], "source": [ "geo_model.set_topography(\n", " source='gdal', filepath=file_path + 'raster7.tif')" ] }, { "cell_type": "markdown", "id": "extreme-drink", "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": "conservative-philip", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:31.521012Z", "start_time": "2022-04-06T09:52:31.187916Z" } }, "outputs": [], "source": [ "gp.plot_2d(geo_model, direction='z', show_lith=False, show_boundaries=False)\n", "plt.grid()" ] }, { "cell_type": "markdown", "id": "db91e4fd", "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": "vocational-intranet", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:31.950110Z", "start_time": "2022-04-06T09:52:31.523010Z" } }, "outputs": [], "source": [ "gp.plot_3d(geo_model, image=False, plotter_type='basic', notebook=True)" ] }, { "cell_type": "markdown", "id": "geological-triple", "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": "external-clerk", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:52:34.750216Z", "start_time": "2022-04-06T09:52:31.953096Z" } }, "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": "above-lending", "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": "impressed-adams", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:54:59.788645Z", "start_time": "2022-04-06T09:52:34.752217Z" } }, "outputs": [], "source": [ "sol = gp.compute_model(geo_model, compute_mesh=True)" ] }, { "cell_type": "code", "execution_count": null, "id": "72de08ed", "metadata": {}, "outputs": [], "source": [ "sol" ] }, { "cell_type": "code", "execution_count": null, "id": "9884651b", "metadata": {}, "outputs": [], "source": [ "geo_model.solutions" ] }, { "cell_type": "markdown", "id": "developing-wallet", "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 layers of the model in the different directions. \n", "\n", "The first section to be plotted is the custom section `Section1` followed by an array of cross sections." ] }, { "cell_type": "code", "execution_count": null, "id": "champion-collectible", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:55:00.758480Z", "start_time": "2022-04-06T09:54:59.790634Z" } }, "outputs": [], "source": [ "gp.plot_2d(geo_model, direction=['x', 'x', 'y', 'y'], cell_number=[25, 75, 25, 75], show_topography=True, show_data=False)" ] }, { "cell_type": "markdown", "id": "6c397cd3", "metadata": {}, "source": [ "Next to the lithology data, we can also plot the calculated scalar field." ] }, { "cell_type": "code", "execution_count": null, "id": "6ff7e5a0", "metadata": {}, "outputs": [], "source": [ "gp.plot_2d(geo_model, direction='y', show_data=False, show_scalar=True, show_lith=False)" ] }, { "cell_type": "markdown", "id": "intimate-connectivity", "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": "developing-comparison", "metadata": { "ExecuteTime": { "end_time": "2022-04-06T09:55:04.041397Z", "start_time": "2022-04-06T09:55:00.760464Z" } }, "outputs": [], "source": [ "gpv = gp.plot_3d(geo_model, image=False, show_topography=True,\n", " plotter_type='basic', notebook=True, show_lith=True)" ] }, { "cell_type": "markdown", "id": "1890840e", "metadata": {}, "source": [ "\n", "\n", "## Post Processing\n", "\n", "Now that the model has been created, we would like to utilize it. In this notebook, you will learn how to sample a raster (interfaces/orientations), how to calculate raster properties and how to perform raster operations and how to visualize input data in **PyVista**.\n", "\n", "\n", "\n", "### Sampling Interfaces and Orientations from Raster\n", "\n", "In the previous tutorials, you have automatically sampled the height information from rasters. In this section you will learn how to use the functions behind the sampling of Z-values to make it more flexible for other purposes. We hereby use the functions `gg.raster.sample_from_array(..)`, `gg.raster.sample_from_raster(..)`, `gg.raster.sample_randomly(..)`. \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/06_sampling_from_rasters.html).\n" ] }, { "cell_type": "code", "execution_count": null, "id": "9c28c05d", "metadata": {}, "outputs": [], "source": [ "extent = [topo_raster.bounds[0], topo_raster.bounds[2], topo_raster.bounds[1], topo_raster.bounds[3]]\n", "extent" ] }, { "cell_type": "markdown", "id": "327e23d0", "metadata": {}, "source": [ "Sampling a single sample. " ] }, { "cell_type": "code", "execution_count": null, "id": "566e1a8d", "metadata": {}, "outputs": [], "source": [ "sample = gg.raster.sample_from_array(array = topo_raster.read(1),\n", " extent = extent,\n", " point_x = 500,\n", " point_y = 500)\n", "\n", "sample" ] }, { "cell_type": "markdown", "id": "e80cf809", "metadata": {}, "source": [ "Sampling from given lists of x and y coordinates." ] }, { "cell_type": "code", "execution_count": null, "id": "b412472e", "metadata": {}, "outputs": [], "source": [ "point_x = [100, 200, 300, 400 ,500]\n", "point_y = [100, 200, 300, 400 ,500]\n", "\n", "sample = gg.raster.sample_from_array(array = topo_raster.read(1),\n", " extent = extent,\n", " point_x = point_x,\n", " point_y = point_y)\n", "\n", "sample" ] }, { "cell_type": "markdown", "id": "b3ccc81e", "metadata": {}, "source": [ "Sampling from given arrays of x and y coordinates." ] }, { "cell_type": "code", "execution_count": null, "id": "59c5e3b3", "metadata": {}, "outputs": [], "source": [ "point_x = np.array([100, 200, 300, 400 ,500.])\n", "point_y = np.array([100, 200, 300, 400 ,500.])\n", "\n", "sample = gg.raster.sample_from_array(array = topo_raster.read(1),\n", " extent = [0, 972, 0, 1069],\n", " point_x = point_x,\n", " point_y = point_y)\n", "\n", "sample" ] }, { "cell_type": "markdown", "id": "5fcba951", "metadata": {}, "source": [ "Sampling a given point from a rasterio object." ] }, { "cell_type": "code", "execution_count": null, "id": "94714c78", "metadata": {}, "outputs": [], "source": [ "sample = gg.raster.sample_from_rasterio(raster = topo_raster,\n", " point_x = 500,\n", " point_y = 500)\n", "\n", "sample" ] }, { "cell_type": "markdown", "id": "c9d52880", "metadata": {}, "source": [ "Sampling from given lists of x and y coordinates." ] }, { "cell_type": "code", "execution_count": null, "id": "f1adcc4f", "metadata": {}, "outputs": [], "source": [ "point_x = [100, 200, 300, 400 ,500]\n", "point_y = [100, 200, 300, 400 ,500]\n", "\n", "sample = gg.raster.sample_from_rasterio(raster = topo_raster,\n", " point_x = point_x,\n", " point_y = point_y)\n", "\n", "sample" ] }, { "cell_type": "markdown", "id": "eb3bcb4b", "metadata": {}, "source": [ "Sampling from given arrays of x and y coordinates." ] }, { "cell_type": "code", "execution_count": null, "id": "5cbc7b8b", "metadata": {}, "outputs": [], "source": [ "point_x = np.array([100, 200, 300, 400 ,500.])\n", "point_y = np.array([100, 200, 300, 400 ,500.])\n", "\n", "sample = gg.raster.sample_from_rasterio(raster = topo_raster,\n", " point_x = point_x,\n", " point_y = point_y)\n", "\n", "sample" ] }, { "cell_type": "markdown", "id": "e9035671", "metadata": {}, "source": [ "Sampling randomly from the raster. This will also return the coordinates of the sampled point as the second entry of the tuple. The same can be done by providing an array and an extent to the function. " ] }, { "cell_type": "code", "execution_count": null, "id": "616e7243", "metadata": {}, "outputs": [], "source": [ "sample = gg.raster.sample_randomly(raster=topo_raster,\n", " n=1)\n", "\n", "sample" ] }, { "cell_type": "code", "execution_count": null, "id": "ac2f1c45", "metadata": {}, "outputs": [], "source": [ "sample = gg.raster.sample_randomly(raster=topo_raster,\n", " n=5)\n", "\n", "sample" ] }, { "cell_type": "markdown", "id": "80b76f28", "metadata": {}, "source": [ "This functionality can be expanded to return a `GeoDataFrame` including the `X`, `Y`, `Z` coordinates, a provided `formation` name and a `crs`. Here, we sample again a `rasterio` object but you can also provide an array and extent to get the same results. By providing a number of `random_samples` instead of x and y coordinates, random samples will be drawn and a `GeoDataFrame` will be returned. \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/08_sampling_interfaces_orientations_from_rasters.html)." ] }, { "cell_type": "code", "execution_count": null, "id": "313255ac", "metadata": {}, "outputs": [], "source": [ "point_x = [100, 200, 300, 400, 500]\n", "point_y = [100, 200, 300, 400, 500]\n", "\n", "gdf = gg.raster.sample_interfaces(raster=topo_raster,\n", " point_x=point_x,\n", " point_y=point_y,\n", " formation='Topography',\n", " crs='EPSG:4326')\n", "\n", "gdf" ] }, { "cell_type": "markdown", "id": "1a3974b9", "metadata": {}, "source": [ "We can not only sample interfaces from a raster but also orientations that returns a `GeoDataFrame` containing also `dip`, `azimuth` and `polarity` values. The same can again be done by sampling from an array or by random sampling. \n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "16343305", "metadata": {}, "outputs": [], "source": [ "point_x = [100, 200, 300, 400 ,500]\n", "point_y = [100, 200, 300, 400 ,500]\n", "\n", "gdf = gg.raster.sample_orientations(raster=topo_raster,\n", " point_x=point_x,\n", " point_y=point_y,\n", " formation='Topography',\n", " crs='EPSG:4326')\n", "\n", "gdf" ] }, { "cell_type": "markdown", "id": "d445f78c", "metadata": {}, "source": [ "\n", "\n", "### Calculating Raster Properties\n", "\n", "In order to calculate the dip and azimuth of orientations, the dip and azimuth of the layer at a respective point needs to be calculated. `GemGIS` has functions implemented to calculate these properties automatically. The slope can be calculated using `gg.raster.calculate_slope(..)`. \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/07_calculating_raster_properties.html)." ] }, { "cell_type": "code", "execution_count": null, "id": "7951491e", "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "\n", "im = plt.imshow(topo_raster.read(1), cmap='gist_earth')\n", "plt.grid()\n", "plt.colorbar(im)" ] }, { "cell_type": "code", "execution_count": null, "id": "b6c033d1", "metadata": {}, "outputs": [], "source": [ "slope = gg.raster.calculate_slope(topo_raster)\n", "im = plt.imshow(slope, cmap='RdYlGn_r', vmin=0, vmax=90)\n", "plt.grid()\n", "cbar = plt.colorbar(im)\n", "cbar.set_label('Slope in degrees')" ] }, { "cell_type": "code", "execution_count": null, "id": "91c54b6b", "metadata": {}, "outputs": [], "source": [ "aspect = gg.raster.calculate_aspect(topo_raster)\n", "im = plt.imshow(aspect, cmap='twilight_shifted')\n", "plt.grid()\n", "cbar = plt.colorbar(im)\n", "cbar.set_label('Gradient in degrees')" ] }, { "cell_type": "markdown", "id": "a548b150", "metadata": {}, "source": [ "\n", "\n", "### Performing Raster Operations\n", "\n", "Different raster operations are implemented in `GemGIS`. These include calculating differences between rasters, resizing rasters, saving rasters and merging rasters. A more detailed tutorial of the functionality can be found [here](https://gemgis.readthedocs.io/en/latest/getting_started/tutorial/09_raster_operations_gemgis.html).\n", "\n" ] }, { "cell_type": "markdown", "id": "8b44269f", "metadata": {}, "source": [ "\n", "\n", "### Visualizing Input Data with PyVista\n", "\n", "The Input Data used for the structural modeling in `GemPy` cannot only be displayed in 2D using `matplotlib` but also in 3D using `PyVista` and implemented functions in `GemGIS`. We use the functions `gg.visualization.create_lines_3d_polydata()` to create PolyData datasets from contour lines, `gg.visualization.create_dem_3d()` to create a StructuredGrid of the digital elevation model, and `gg.visualization.create_points_3d()` to PolyData datasets of interface points. \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/10_visualizing_data_with_pyvista.html)." ] }, { "cell_type": "code", "execution_count": null, "id": "33eaff50", "metadata": {}, "outputs": [], "source": [ "lines = gg.visualization.create_lines_3d_polydata(gdf=topo)\n", "lines" ] }, { "cell_type": "code", "execution_count": null, "id": "4d3c2c81", "metadata": {}, "outputs": [], "source": [ "grid = gg.visualization.create_dem_3d(dem=topo_raster)\n", "grid" ] }, { "cell_type": "code", "execution_count": null, "id": "29ee74c6", "metadata": {}, "outputs": [], "source": [ "points_mesh = gg.visualization.create_points_3d(gdf=interfaces_coords)\n", "points_mesh" ] }, { "cell_type": "code", "execution_count": null, "id": "0ee24638", "metadata": {}, "outputs": [], "source": [ "sargs = dict(fmt=\"%.0f\", color='black')\n", "\n", "p = pv.Plotter(notebook=True)\n", "\n", "p.add_mesh(mesh=lines, color='black')\n", "p.add_mesh(grid, scalar_bar_args=sargs, clim=[500,1500])\n", "p.add_mesh(mesh=points_mesh, color='red')\n", "\n", "p.show_grid(color='black')\n", "p.set_background(color='white')\n", "p.show()" ] }, { "cell_type": "markdown", "id": "f1667ce8", "metadata": {}, "source": [ "\n", "# Conclusions\n", "\n", "
\n", "In this tutorial, you have learnt the following:
\n", "- How to build your fifth GemPy model with input data generated through GemGIS
\n", "- How to sample rasters (interfaces/orientations)
\n", "- How to calculate raster properties and how to perform raster operations
\n", "- How to visualize the input data in PyVista
\n", "
\n", "\n", "\n", "# Outlook\n", "\n", "
\n", "In the next tutorial, you will learn the following:
\n", "- How to build your sixth GemPy model with input data generated through GemGIS
\n", "- How to create orientations from isolines on map
\n", "\n", "\n", "
\n", "\n", "[Take me to the next notebook on Github](https://nbviewer.org/github/cgre-aachen/gemgis_data/blob/main/notebooks/03_folded_layers/example03_folded_layers.ipynb)\n", "\n", "[Take me to the next notebook locally](example03_folded_layers.ipynb)\n", "\n", "\n", "Source: Powell, D. (1995): Interpretation geologischer Strukturen durch Karten - Eine praktische Anleitung mit Aufgaben und Lösungen, page 31, figure 26 A, Springer Verlag Berlin, Heidelberg, New York, ISBN: 978-3-540-58607-4.\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 }