Welcome to the course "From Maps to Models - Tutorials for structural geological modeling using GemPy and GemGIS". In this course, you will learn how to obtain and process spatial data with the Python open-source package GemGIS, how to create structural geological models using the Python open-source package GemPy and how to post process the models to obtain certain information from your models. Both packages were developed and are maintained by departments of Prof. Florian Wellmann at RWTH Aachen University and Fraunhofer IEG. The course material is adapted from two undergraduate courses at RWTH Aachen University, Germany (Geological methods and geological map & Introduction to geological maps) taught by Prof. Stefan Back and Prof. Susanne Buiter. The tutorials are part of the publication by Jüstel et al. (2023) in The Journal of Open Source Education.
The course is divided into four sections. The main part is further divided into 6 sections
In this notebook, the Introduction to the course
and Introduction to structural geological modeling using GemPy and GemGIS
are covered.
The tutorials provided in this repository aim at enabling students who want to become familiar with structural geological modeling in GemPy
to have a smooth start. Further, the tutorials serve as aids for researchers to set up and parametrize their structural geological models properly. Lectures are encouraged to use and adapt the provided tutorials for their own Python and Structural Geological Modeling classes. In addition, the pre- and post-processing tasks introduced in several notebooks outline the possibilities of GemGIS
apart from using it in conjunction with GemPy
.
The students should have basic knowledge in creating vector data (shape files) in geoinformation systems like QGIS and basic knowledge of programming in Python using Jupyter Notebooks.
Source: https://gisanalyse.de/category/gis-analysesThe following learning activities are recommended during the course specified for each week for a 14 weeks period of lectures and exercises. The students may choose more than the three required models by themselves if more than 3 models are available for each section.
The easiest way to access the modeling notebooks and data is to download the repository including all data from Github. The following folders are then needed:
If you have questions about the tutorials or the modeling with GemPy
or data processing with GemGIS
, please open a new Discussion in GemPy Discussions or GemGIS Discussions
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
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) and Bennison (1990). References for maps with unknown origin will gladly be added.
The following sections provide an introduction to geological maps, geological models, the Python open-source packages GemPy
und GemGIS
, installation instructions to create a running conda
environment and further resources to get started with GemPy
and GemGIS
.
Geological maps have been used since egyptian times to illustrate the spatial distribution of rocks on the surface and in the shallow subsurface. Geological maps are used to illustrate the distribution of resources (minerals, building materials), to outline geohazards, geoheritages and many more purposes. Geological maps can be created on various scales (e.g., 1:10,000, 1:100,000) or from various scales (work directly in the field or through remote sensing techniques). Additional information from drillings, geochemical investigations or geophysical investigations can be integrated to complement the maps. Often, the resulting map shows a combination of the age of the rocks (oldest at the bottom), the stratigraphic column, and the lithology of the rocks (what kind of rock it is) as well as the main tectonic features of the respective region.
Depending on the location on earth or on extraterrestrial objects (Moon, Mars), geological maps represent projections of the spherical surface of the object in 3D onto a 2D map in a specific coordinate range. Different projections exist but their introducions are beyond the scope of this course now as you will work with small scale local examples where projections are irrelevant. Nevertheless, it is important to mention that projections are always distorting the true features of the spherical surface when projecting them to 2D. A common projection is the cylindrical Universal-Transverse-Mercator (UTM) projection.
Apart from arid regions like Svalbard (Spitzbergen) or deserts, the geology is mostly covered by erosional deposits (alluvial, fluvial, aeolian) which in itself make up a new geological map. Further vegetation or civilization (cities) may prevent us from seeing the outcropping rocks. A geologist can therefore not always comprehent the full picture of the subsurface structures and needs to interpolate between the available information and update his/her map when new information becomes available (new outcrop, drilling, construction work). Therefore, different geologists may have different interpretations which may all be true given the available data and following basic geological principles as introduced below.
Source: Lecture material of Susanne Buiter @ RWTH Aachen University, Waldron, J.; Snyder, J. (2020): Geological Structures: a practical introduction. University of Alberta, https://doi.org/10.29173/oer3Geological information can be visualized in different ways. Geological maps (view from the top considering the topography) is one way. A second way is to show the geological structures in a 3D structural geological model which you are going to create in this tutorial series and work with. Two more ways of visulazing geological data are using cross sections and depth maps which are mostly a result of the geological map or a 3D structural model.
In this course you will also learn how to read geological expressions at the surface and what they indicate for the course of the layers in the subsurface. What is the orientation of the dotted layer in the different images below? Horizontal? Vertical? Inclined?
Source: Lecture material of Susanne Buiter @ RWTH Aachen University, Borradaile, G. (2014): Understanding Geology Through Maps, Elsevier, Amsterdam, Netherlands, ISBN: 978-0-12-800866-9.Structural geological models (sensu Wellmann & Caumon, 2018) represent geometric elements in the subsurface, so-called surfaces (Fig. 1). These surfaces may be lithological boundaries such as a boundary between a sandstone or a claystone, which are used mostly in this course, stratigraphic boundaries between layers of different age, which are mostly used in the real world, faults separating a continuous layer or intrusions. A layer is defined here as rock mass exhibiting the same properties such as lithology or stratigraphic age. A layer is bound by a top and a bottom surface. The structural geological model is created in the first place according to the stratigraphic principles of Steno (1669) that more recent strata are deposited roughly horizontally on top of older strata and that variations such as folded, faulted or unconformable layers are caused by geological event.
Fig. 1: Geological model of the Perth basin (Australia) rendered using GemPy
on the in-built Python in Blender, spheres and cones represent the input data. From de la Varga et al. (2019).
Source: Powell, D. (1995): Interpretation geologischer Strukturen durch Karten - Eine praktische Anleitung mit Aufgaben und Lösungen, page 14, figure 7, Springer Verlag Berlin, Heidelberg, New York, ISBN: 978-3-540-58607-4.
The course of the layers can be observed on the surface in outcrops or quarries (Fig. 2). Further, direct observations of layer boundaries can be made in the subsurface in boreholes or indirectly using geophysical methods such as seismic representing the seismic impedance (product of rock density and seismic velocity). This course focuses on observations on the surface and in boreholes.
Fig. 2: Quarry Binsfeldhammer, Stolberg, Aachen City Region, Germany.
The information that can be obtained from observed layer boundaries are the locations of these boundaries (x
,y
,z
) and which boundary they represent (for this course termed formation
), and the orientation of this layer at a measured location expressed as dip
and azimuth
(Fig. 3). The locations should be expressed in meters (m) rather than geographical coordinates (longitude/latitude). The dip and azimuth are expressed 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).
Fig. 3: Example of strike and dip on tilted sedimentary beds.
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
Structural geological models describe the spatial distribution of layers in the subsurface. This includes the extent and depth distribution of either lithological layers (e.g. sand or clay, sandstones, carbonates, evaporites, etc.) or stratigraphic units for more regional questions of the structure of the subsurface (e.g. Base Tertiary, transition to basement rocks, etc.). The resulting structures can be used to parametrize the rocks between these constructed bounding surfaces for a wide range of applications. These include but are not limited to hydrogeological applications, geohazards, ressource extraction (e.g. fresh water, thermal waters, hydrocarbons, minerals), storage applications (e.g. thermal waters, gases such as CO2 or hydrogen, nuclear waste).
GemPy
and GemGIS
are two Python open-source packages developed and maintained by departments of Prof. Florian Wellmann at RWTH Aachen University and Fraunhofer IEG. GemPy
is capable of creating structural geological models based on an implicit modeling approach using universal co-kriging (de la Varga et al., 2019) and the potential field approach by Lajaunie et al. (1997).
A simple synclinal model is shown to illustrate the principles of GemPy
(Fig. 4, left). The first input data used for the interpolation of geological surfaces, here stratigraphic surfaces, are the location (x, y, z) of stratigraphic boundaries mapped in surface outcrops, interpolated on geological maps, in wells or seismic data in the subsurface, or constrained from additional geophysical methods (e.g. gravity, ERT). Orientation data (x, y, z, dip, azimuth) obtained from outcrop measurements or reference values from literature are the second type of input data (Fig. 4, left. colored arrows). Subparallel layers/formations are combined in one series and are interpolated using one scalar field/potential field where a certain isovalue of this field corresponds to the location of a stratigraphic boundary (Fig. 4, right). The measured orientations represent the gradient of the field (Fig. 4, right). Layers that are separated by an unconformity are modeled in separate series. The same accounts for faults where each fault is modeled as its own series. The order of the series and the formations from old at the bottom and young at the top with faults above the series and formations define the stratigraphic pile. This will be introduced in more detail during the modeling of basic structures.
The interpolation in GemPy
is meshless working with a regular grid where the resulting layer or stratigraphic unit is assigned to the respective grid cell (voxel). The actual resulting structural geological model is therefore a so-called lith_block
represented by a NumPy array. A marching-cube algorithm is applied to calculate the 2.5D surface meshes that can be visualized in PyVista (Sullivan and Kaszynski, 2019) and that can be used for post-processing in GemGIS
.
Fig. 4: Structural geological model showing the input data points, input orientations (arrows) and resulting meshes (continuous colored lines) with the resulting lith_block
(left) and scalar field (right).
The GemGIS
package has been developed as an Add-On to GemPy
to accelerate and simplify the processing of spatial data as input data for the structural geological modeling with GemPy
(Fig 5). It is wrapping and extends the functionality of already well-known spatial data open-source Python packages such as GeoPandas (Bosche et al., 2022), Rasterio (Gillies et al., 2019) or PyVista. In addition, it provides functions to utilize the created structural geological models further (creating virtual boreholes, creating depth maps, etc.).
Fig. 5: From Maps to Models. Source: Powell, D. (1995): Interpretation geologischer Strukturen durch Karten - Eine praktische Anleitung mit Aufgaben und Lösungen, page 15, figure 10 A, Springer Verlag Berlin, Heidelberg, New York, ISBN: 978-3-540-58607-4.
The following resources are recommended for GemPy
and GemGIS
. In addition, it is recommended to join the Software Underground (SWUNG) Slack Workspace. This is "the place for scientists and engineers that love rocks and computers". In this community with more than 4000 enthusiasts, there are dedicated channels for the structural geological modeling with #gempy, #geospatial data, #open-geosciences and many more. The annual "Transform" conference is the official conference of the Software Underground with many tutorial sessions introducing the different open-source packages that have emerged from this community among other presentations and lightning talks.
For extensive installation instructions, see GemPy Installation and GemGIS Installation.
Both packages are available on PyPi and conda:
In this section, we will provide a short version of the installation instructions to get you started.
It is recommended to use the latest Anaconda distribution.
It is recommended to create a new virtual environment when using GemGIS
to avoid conflicts with already existing projects. This step and all following steps will be performed within the Anaconda Prompt (Anaconda3).
Creating a new environment in Anaconda with fixed Python version:
conda create -n gemgis python==3.10
Activate the virtual environment:
conda activate gemgis
The gemgis environment now replaced the base environment which is indicated by gemgis
in front of your path.
Several packages need to be installed in order to use GemGIS
. It is recommended to use Jupyter Notebooks
when working with GemGIS
. Packages are installed using the conda-forge
channel.
Install Jupyter Notebooks:
conda install -c conda-forge jupyter
You can start a new kernel by executing jupyter notebook
in your Anaconda Prompt. A new kernel will then open in your browser.
GemGIS
and all its dependencies can be installed via conda-forge:
conda install -c conda-forge gemgis
Two of the main packages that GemGIS
is dependent on are rasterio and GeoPandas. It is recommended to install these packages separately as they both depend on the GDAL translator library for raster and vector geospatial data. In addition, many smaller libraries like shaply or fiona will also be installed properly.
Install the latest versions of GeoPandas and Rasterio (as of 2023-01-01). Please mind the quotation marks that are necessary when specifying the version numbers:
conda install -c conda-forge geopandas">=0.12.2" rasterio">=1.3.4" pyvista">=0.38.2"
pip install gemgis
Once GemGIS
has successfully been installed, GemPy
can be installed via pip. Please mind that GemPy
is an optional dependency of GemGIS
.
pip install gempy
If you want to use GemPy
and/or GemGIS
, do not forget to make the correct contributions.
For GemPy, please cite the publication by de la Varga et al. (2019):
@article{delaVarga2019,
author = {de la Varga, M. and Schaaf, A. and Wellmann, F.},
title = {GemPy 1.0: open-source stochastic geological modeling and inversion},
journal = {Geoscientific Model Development},
volume = {12},
year = {2019},
number = {1},
pages = {1-32},
url = {https://gmd.copernicus.org/articles/12/1/2019/},
doi = {https://doi.org/10.5194/gmd-12-1-2019}
}
For GemGIS, please cite the publication by Jüstel et al. (2022):
@article{Jüstel2022,
doi = {https://doi.org/10.21105/joss.03709},
url = {https://doi.org/10.21105/joss.03709},
year = {2022},
publisher = {The Open Journal},
volume = {7},
number = {73},
pages = {3709},
author = {Jüstel, A. and Endlein Correira, A. and Pischke, P. and de la Varga, M. and Wellmann, F.},
title = {GemGIS - Spatial Data Processing for Geomodeling},
journal = {Journal of Open Source Software}
}
In the next four notebooks, the principles of GemPy
will be explained. Here, you will model planar layers, folded layers, faulted layers and unconformable layers truncated by an unconformity. No topography is considered in these models.
Bennison, G. M. (1998): An Introduction to Geological Structures and Maps. Boston, MA: Springer US. https://link.springer.com/book/10.1007/978-1-4615-9630-1
Borradaile, G. (2014): Understanding Geology Through Maps, Elsevier, Amsterdam, Netherlands, ISBN: 978-0-12-800866-9
Bossche, M. F.; McB., James; Wasserman, J.; Richards, M; et al. (2022): geopandas/geopandas: v0.12.2: Zenodo. https://doi.org/10.5281/zenodo.2585848
de La Varga, Miguel ; Schaaf, Alexander; Wellmann, Florian (2019): GemPy 1.0: open-source stochastic geological modeling and inversion. 0001_delaVarga2019_GemPy.pdf. In: Geosci. Model Dev. 12 (1), S. 1–32. https://doi.org/10.5194/gmd-12-1-2019.
Gillies, S.; et al. (2019): Rasterio: geospatial raster I/O for Python programmers. https://github.com/rasterio/rasterio.
Jüstel, A.; Endlein Correira, A.; Pischke, M.; de La Varga, M.; Wellmann, F. (2022): GemGIS - Spatial Data Processing for Geomodeling. In: JOSS 7 (73), S. 3709. https://doi.org/10.21105/joss.03709.
Lajaunie, Christian; Courrioux, Gabriel; Manuel, Laurent (1997): Foliation fields and 3D cartography in geology: Principles of a method based on potential interpolation. In: Math Geol 29 (4), S. 571–584. https://doi.org/10.1007/BF02775087.
Powell, Derek (1995): Interpretation geologischer Strukturen durch Karten. Eine praktische Anleitung mit Aufgaben und Lösungen. Berlin, Heidelberg, New York, Barcelona, Budapest, Hong Kong, London, Mailand, Paris, Tokyo: Springer. https://link.springer.com/book/9783540586074.
Steno, N. (1669): De solido intra solidium naturaliter contento dissertationis prodromus. Florence.
Sullivan, C.; Kaszynski, Alexander (2019): PyVista: 3D plotting and mesh analysis through a streamlined interface for the Visualization Toolkit (VTK). In: JOSS 4 (37), S. 1450. DOI: https://doi.org/10.21105/joss.01450.
Waldron, J.; Snyder, J. (2020): Geological Structures: a practical introduction. University of Alberta, https://doi.org/10.29173/oer3
Wellmann, Florian; Caumon, Guillaume (2018): 3-D Structural geological models: Concepts, methods, and uncertainties. In: Cedric Schmelzbach (Hg.): Advances in Geophysics, Bd. 59: Elsevier (Advances in Geophysics), S. 1–121. https://doi.org/10.1016/bs.agph.2018.09.001.
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
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) and Bennison (1990). References for maps with unknown origin will gladly be added.