Merge pull request #255 from cgre-aachen/dev_gemgis2

GemGIS 1.0.8

docs/conf.py

@@ -24,7 +24,7 @@
 author = 'Alexander Juestel'
 
 # The full version, including alpha/beta/rc tags
-release = '1.0.7'
+release = '1.0.8'
 version = release
 
 # -- GemGIS configuration ---------------------------------------------------

gemgis/__init__.py

@@ -20,11 +20,11 @@
 
 __version_date__ = '2023-02-07'
 
-__version__ = '1.0.7'
+__version__ = '1.0.8'
+
+__changelog__ = """What is new in version 1.0.8:
 
-__changelog__ = """What is new in version 1.0.7:
 
-- Fix for #249
 
 """
 

gemgis/postprocessing.py

@@ -24,16 +24,36 @@
 import numpy as np
 from typing import List, Union
 from gemgis import gemgis
+import pyvista as pv
+import pyproj
 
 
-def extract_lithologies(geo_model, extent, crs):
+def extract_lithologies(geo_model,
+                        extent: list,
+                        crs: Union[str, pyproj.crs.crs.CRS]) -> gpd.geodataframe.GeoDataFrame:
+    """Extracting the geological map as GeoDataFrame
+
+    Parameters:
+    ___________
+
+        geo_model: gp.core.model.Project
+            GemPy geo_model
+
+        extent: list
+            Extent of geo_model
+
+        crs: Union[str, pyproj.crs.crs.CRS]
+            Coordinate References System
+
+    """
     # Trying to import matplotlib but returning error if matplotlib is not installed
     try:
         import matplotlib.pyplot as plt
     except ModuleNotFoundError:
         raise ModuleNotFoundError(
             'Matplotlib package is not installed. Use pip install matplotlib to install the latest version')
 
+    # Trying to import gempy but returning error if gempy is not installed
     try:
         import gempy as gp
     except ModuleNotFoundError:
@@ -84,8 +104,9 @@ def extract_lithologies(geo_model, extent, crs):
             fm.append(fm_name)
             geo.append(poly)
 
-    lith = gpd.GeoDataFrame({"formation": fm}, geometry=geo)
-    lith.crs = crs
+    lith = gpd.GeoDataFrame({"formation": fm},
+                            geometry=geo,
+                            crs=crs)
 
     return lith
 
@@ -319,4 +340,167 @@ def save_model(geo_model, path):
     np.save(path + "02_extent.npy", geo_model.grid.regular_grid.extent)
     np.save(path + "03_resolution.npy", geo_model.grid.regular_grid.resolution)
 
+
+def extract_orientations_from_mesh(mesh: pv.core.pointset.PolyData,
+                                   crs: Union[str, pyproj.crs.crs.CRS]) -> gpd.geodataframe.GeoDataFrame:
+    """Extracting orientations (dip and azimuth) from PyVista Mesh
+
+    Parameters:
+    ___________
+
+        mesh: pv.core.pointset.PolyData
+            PyVista Mesh from which the orientations will be extracted
+
+        crs: Union[str, pyproj.crs.crs.CRS]
+            Coordinate reference system of the returned GeoDataFrame
+
+    Returns:
+    ________
+
+        gdf_orientations: gpd.geodataframe.GeoDataFrame
+            GeoDataFrame consisting of the
+
+    """
+
+    # Checking that the provided mesh is of type Polydata
+    if not isinstance(mesh, pv.core.pointset.PolyData):
+        raise TypeError('Mesh must be provided as PyVista Polydata')
+
+    # Checking that the provided mesh if of type string or a pyproj CRS object
+    if not isinstance(crs, (str, pyproj.crs.crs.CRS)):
+        raise TypeError('CRS must be provided as string or pyproj CRS object')
+
+    # Computing the normals of the mesh
+    mesh_normals = mesh.compute_normals()
+
+    # Calculating the dips
+    dips = [90 - np.rad2deg(-np.arcsin(mesh_normals['Normals'][i][2])) * (-1) for i in
+            range(len(mesh_normals['Normals']))]
+
+    # Calculating the azimuths
+    azimuths = [np.rad2deg(np.arctan(mesh_normals['Normals'][i][0] / mesh_normals['Normals'][i][1])) + 180 for i in
+                range(len(mesh_normals['Normals']))]
+
+    # Getting cell centers
+    points_z = [geometry.Point(point) for point in mesh.cell_centers().points]
+
+    # Creating GeoDataFrame
+    gdf_orientations = gpd.GeoDataFrame(geometry=points_z,
+                                        crs=crs)
+
+    # Appending X, Y, Z Locations
+    gdf_orientations['X'] = mesh.cell_centers().points[:, 0]
+    gdf_orientations['Y'] = mesh.cell_centers().points[:, 1]
+    gdf_orientations['Z'] = mesh.cell_centers().points[:, 2]
+
+    # Appending dips and azimuths
+    gdf_orientations['dip'] = dips
+    gdf_orientations['azimuth'] = azimuths
+
+    return gdf_orientations
+
+
+def calculate_dip_and_azimuth_from_mesh(mesh: pv.core.pointset.PolyData) -> pv.core.pointset.PolyData:
+    """Calculating dip and azimuth values for a mesh and setting them as scalars for subsequent plotting
+
+    Parameters:
+    ___________
+
+        mesh: pv.core.pointset.PolyData
+            PyVista Mesh for which the dip and the azimuth will be calculated
+
+    Returns:
+    ________
+
+        mesh: pv.core.pointset.PolyData
+            PyVista Mesh with appended dips and azimuths
+
+    """
+
+    # Checking that the provided mesh is of type Polydata
+    if not isinstance(mesh, pv.core.pointset.PolyData):
+        raise TypeError('Mesh must be provided as PyVista Polydata')
+
+    # Computing the normals of the mesh
+    mesh.compute_normals(inplace=True)
+
+    # Calculating the dips
+    dips = [90 - np.rad2deg(-np.arcsin(mesh['Normals'][i][2])) * (-1) for i in
+            range(len(mesh['Normals']))]
+
+    # Calculating the azimuths
+    azimuths = [np.rad2deg(np.arctan(mesh['Normals'][i][0] / mesh['Normals'][i][1])) + 180 for i in
+                range(len(mesh['Normals']))]
+
+    # Assigning dips and azimuths to scalars
+    mesh['Dips [°]'] = dips
+    mesh['Azimuths [°]'] = azimuths
+
+    return mesh
+
+
+def crop_block_to_topography(geo_model) -> pv.core.pointset.UnstructuredGrid:
+    """Cropping GemPy solutions block to topography
+
+    Parameters:
+    ___________
+
+        geo_model: gp.core.model.Project
+
+    Returns:
+    ________
+
+        grid: pv.core.pointset.UnstructuredGrid
+
+
+    """
+
+    # Trying to import GemPy
+    try:
+        import gempy as gp
+    except ModuleNotFoundError:
+        raise ModuleNotFoundError(
+            'GemPy package is not installed. Use pip install gempy to install the latest version')
+
+    # Trying to import PVGeo
+    try:
+        from PVGeo.grids import ExtractTopography
+    except ModuleNotFoundError:
+        raise ModuleNotFoundError(
+            'PVGeo package is not installed. Use pip install pvgeo to install the lastest version')
+
+    # Creating StructuredGrid
+    grid = pv.UniformGrid()
+
+    # Setting Grid Dimensions
+    grid.dimensions = np.array(geo_model.solutions.lith_block.reshape(geo_model.grid.regular_grid.resolution).shape) + 1
+
+    # Setting Grid Origin
+    grid.origin = (geo_model.grid.regular_grid.extent[0],
+                   geo_model.grid.regular_grid.extent[2],
+                   geo_model.grid.regular_grid.extent[4])
+
+    # Setting Grid Spacing
+    grid.spacing = ((geo_model.grid.regular_grid.extent[1] - geo_model.grid.regular_grid.extent[0]) /
+                    geo_model.grid.regular_grid.resolution[0],
+                    (geo_model.grid.regular_grid.extent[3] - geo_model.grid.regular_grid.extent[2]) /
+                    geo_model.grid.regular_grid.resolution[1],
+                    (geo_model.grid.regular_grid.extent[5] - geo_model.grid.regular_grid.extent[4]) /
+                    geo_model.grid.regular_grid.resolution[2])
+
+    # Setting Cell Data
+    grid.cell_data['values'] = geo_model.solutions.lith_block.reshape(geo_model.grid.regular_grid.resolution).flatten(
+        order='F')
+
+    # Creating Polydata Dataset
+    topo = pv.PolyData(geo_model._grid.topography.values)
+
+    # Interpolating topography
+    topo.delaunay_2d(inplace=True)
+
+    extracted = ExtractTopography(tolerance=5,
+                                  remove=True).apply(grid, topo)
+
+    return extracted
+
 # TODO: Create function to export qml layer from surface_color_dict

gemgis/vector.py

@@ -2153,9 +2153,7 @@ def explode_linestring_to_elements(linestring: shapely.geometry.linestring.LineS
         raise ValueError('LineString must contain at least two vertices')
 
     # Splitting the LineString into single elements and returning a list of LineStrings
-    splitted_linestrings = [list(ops.split(list(ops.split(linestring, geometry.Point(linestring.coords[i + 1])).geoms)[0],
-                                      geometry.Point(linestring.coords[i])).geoms)[-1]
-                            for i in range(len(linestring.coords) - 1)]
+    splitted_linestrings = list(map(shapely.geometry.linestring.LineString, zip(linestring.coords[:-1], linestring.coords[1:])))
 
     return splitted_linestrings
 

gemgis/visualization.py

@@ -34,6 +34,7 @@
 from shapely.geometry import LineString
 import pyproj
 import matplotlib
+import matplotlib.colors as mcolors
 import matplotlib.pyplot as plt
 
 
@@ -116,7 +117,7 @@ def create_lines_3d_polydata(gdf: gpd.geodataframe.GeoDataFrame) -> pv.core.poin
 
     # Creating PyVista PolyData containing the lines and vertices
     poly = pv.PolyData()
-    poly.lines = lines
+    poly.lines = np.array(lines)
     poly.points = points
 
     return poly
@@ -4182,3 +4183,112 @@ def get_batlow_cmap() -> matplotlib.colors.ListedColormap:
     ax.set_axis_off()
 
     return cmap_batlow
+
+def get_color_lot(geo_model,
+                  lith_c: pd.DataFrame = None,
+                  index='surface',
+                  is_faults: bool = True,
+                  is_basement: bool = False) -> pd.Series:
+    """Method to get the right color list depending on the type of plot.
+       Borrowed from https://github.com/cgre-aachen/gempy/blob/6aed72a4dfa26830df142a0461294bd9d21a4fa4/gempy/plot/vista.py#L133-L167
+
+    Parameters
+    __________
+
+        geo_model : gp.core.model.Project
+            Previously calculated GemPy Model
+
+        lith_c : pd.DataFrame
+            Pandas series with index surface names and values hex strings with the colors
+
+        index : str
+            Index provided as string
+
+        is_faults : bool
+            Return the colors of the faults. This should be true for surfaces and input data and false for scalar values.
+
+        is_basement : bool
+            Return or not the basement. This should be true for the lith block and false for surfaces and input data.
+    """
+    if lith_c is None:
+        surf_df = geo_model._surfaces.df.set_index(index)
+        unique_surf_points = np.unique(geo_model._surface_points.df['id']).astype(int)
+
+        if len(unique_surf_points) != 0:
+            bool_surf_points = np.zeros(surf_df.shape[0], dtype=bool)
+            bool_surf_points[unique_surf_points.astype('int') - 1] = True
+
+            surf_df['isActive'] = (surf_df['isActive'] | bool_surf_points)
+
+            if is_faults is True and is_basement is True:
+                lith_c = surf_df.groupby('isActive').get_group(True)['color']
+            elif is_faults is True and is_basement is False:
+                lith_c = surf_df.groupby(['isActive', 'isBasement']).get_group((True, False))['color']
+            else:
+                lith_c = surf_df.groupby(['isActive', 'isFault']).get_group((True, False))[
+                    'color']
+
+    color_lot = lith_c
+
+    return color_lot
+
+
+def get_mesh_geological_map(geo_model) -> Tuple[pv.core.pointset.PolyData,
+                                                matplotlib.colors.ListedColormap,
+                                                bool]:
+    """Getting the geological map of a GemPy Model draped over the topography as mesh.
+       Borrowed from https://github.com/cgre-aachen/gempy/blob/6aed72a4dfa26830df142a0461294bd9d21a4fa4/gempy/plot/vista.py#L512-L604
+
+    Parameters:
+    ___________
+
+        geo_model : gp.core.model.Project
+            Previously calculated GemPy Model
+
+    Returns:
+    ________
+
+        polydata: pv.core.PolyData
+            PyVista Mesh containing the geological map draped over the topography
+
+        cm : matplotlib.colors.ListedColormap
+            Colormap for plotting
+
+        rgb : bool
+            Boolean to use rgb=True when plotting
+
+    """
+
+    # Getting topography values from geo_model
+    topography = geo_model._grid.topography.values
+
+    # Creating polydata dataset
+    polydata = pv.PolyData(topography)
+
+    # Getting color values
+    colors_hex = get_color_lot(geo_model=geo_model,
+                               is_faults=False,
+                               is_basement=True, index='id')
+
+    colors_rgb_ = colors_hex.apply(lambda val: list(mcolors.hex2color(val)))
+    colors_rgb = pd.DataFrame(colors_rgb_.to_list(),
+                              index=colors_hex.index) * 255
+
+    sel = np.round(geo_model.solutions.geological_map[0]).astype(int)[0]
+
+    # Converting color values
+    scalars_val = pv.convert_array(colors_rgb.loc[sel].values,
+                                   array_type=3)
+
+    # Creating colormap
+    cm = mcolors.ListedColormap(list(get_color_lot(is_faults=True,
+                                                   is_basement=True)))
+    rgb = True
+
+    # Interpolating the polydata and assigning values
+    polydata.delaunay_2d(inplace=True)
+    polydata['id'] = scalars_val
+    polydata['height'] = topography[:, 2]
+
+    return polydata, cm, rgb
+

setup.py

@@ -1,7 +1,7 @@
 from setuptools import setup, find_packages
 from os import path
 
-version = '1.0.7'
+version = '1.0.8'
 
 # Loading Readme for Description on PyPi
 this_directory = path.abspath(path.dirname(__file__))