Merge pull request #499 from cgre-aachen/development_janN

[DOC] Warning of currently incomaptible pandas version and export functionality for shemat-suite

Author: Leguark

README.md

@@ -11,6 +11,9 @@
 [![DOI](https://zenodo.org/badge/96211155.svg)](https://zenodo.org/badge/latestdoi/96211155)
 [![DOCKER](https://img.shields.io/docker/cloud/automated/leguark/gempy.svg)](https://cloud.docker.com/repository/docker/leguark/gempy)
 
+:warning: **Warning: GemPy requires pandas version < 1.4.0. The new pandas release is not compatible with GemPy.  
+    We're actively working on this issue for a future release.  
+Please make sure to use Pandas version 1.3.x when working with GemPy for the time being.** :warning:
 ## Overview
 
 [GemPy](https://www.gempy.org/) is a Python-based, **open-source geomodeling library**. It is
@@ -66,6 +69,8 @@ Follow these [guidelines](https://github.com/cgre-aachen/gempy/blob/WIP_readme-u
 https://www.sciencedirect.com/science/article/pii/B9780128140482000156).
 In Developments in Structural Geology and Tectonics (Vol. 5, pp. 189-204). Elsevier.
 
+A continuously growing list of gempy-applications (e.g. listing real-world models) can be found [here](https://hackmd.io/@Japhiolite/B1juPvCxc).
+
 ## Gallery
 
 ### Geometries

gempy/__init__.py

@@ -6,6 +6,7 @@
 """
 import sys
 import os
+import pandas
 
 import warnings
 
@@ -37,6 +38,10 @@
 from gempy.plot import _plot as _plot
 
 assert sys.version_info[0] >= 3, "GemPy requires Python 3.X"  # sys.version_info[1] for minor e.g. 6
+assert pandas.__version__ <= '1.4.0', \
+    "GemPy requires pandas version < 1.4.0. The new pandas release is not compatible with GemPy.\n" \
+    "We're actively working on this issue for a future release.\n"
+
 __version__ = '2.2.10'
 
 if __name__ == '__main__':

gempy/assets/topology.py

@@ -53,6 +53,18 @@ def compute_topology(
         n_shift: int = 1,
         voxel_threshold: int = 1
 ):
+    """Compute topology of a GemPy model
+
+    Args:
+        geo_model (Project): GemPy model project
+        cell_number (int, optional): Cell number in chosen direction. Defaults to None.
+        direction (str, optional): one of the model's dimensions, x, y, or z. Defaults to None.
+        n_shift (int, optional): number of voxels the model is shifted and then substracted for finding interfaces. Defaults to 1.
+        voxel_threshold (int, optional): amount of voxels which have to be connected to be considered into the topology calculation. Defaults to 1.
+
+    Returns:
+        edges, centroids [numpy array]: edges and centroids of the topology graph
+    """
     res = geo_model._grid.regular_grid.resolution
     fb = _get_fb(geo_model)
     lb = _get_lb(geo_model)

gempy/core/data_modules/geometric_data.py

@@ -3,6 +3,7 @@
 import warnings
 from typing import Union, Iterable
 
+import numpy
 import numpy as np
 import pandas as pn
 
@@ -291,6 +292,19 @@ def add_surface_points(self, x: Union[float, np.ndarray], y: Union[float, np.nda
 
             self.df.loc[idx, ['X', 'Y', 'Z']] = coord_array.astype('float64')
             self.df.loc[idx, 'surface'] = surface
+            # if isinstance(surface, np.ndarray): # If surface is a numpy array
+            #     # self.df.loc[idx, 'surface'] = surface
+            #     if self.df['surface'].dtype == 'category':
+            #         for s in surface:
+            #             if s not in self.df['surface'].cat.categories:
+            #                 self.df['surface'].cat.add_categories(s, inplace=True)
+            #             self.df.loc[idx, 'surface'] = s
+            #     else:
+            #         for s in surface:
+            #             self.df.loc[idx, 'surface'] = s
+            #
+            # else:
+            #     self.df.loc[idx, 'surface'] = surface
         # ToDO test this
         except ValueError as error:
             self.del_surface_points(idx)
@@ -591,6 +605,11 @@ def add_orientation(self, x, y, z, surface, pole_vector: Union[list, tuple, np.n
 
         if pole_vector is not None:
             self.df.loc[idx, ['X', 'Y', 'Z', 'G_x', 'G_y', 'G_z']] = np.array([x, y, z, *pole_vector], dtype=float)
+            # if type(surface) is numpy.ndarray:
+            #     for s in surface:
+            #         self.df.loc[idx, 'surface'] = s
+            # else:
+            #     self.df.loc[idx, 'surface'] = surface
             self.df.loc[idx, 'surface'] = surface
 
             self.calculate_orientations(idx)
@@ -601,6 +620,11 @@ def add_orientation(self, x, y, z, surface, pole_vector: Union[list, tuple, np.n
             if orientation is not None:
                 self.df.loc[idx, ['X', 'Y', 'Z', ]] = np.array([x, y, z], dtype=float)
                 self.df.loc[idx, ['azimuth', 'dip', 'polarity']] = np.array(orientation, dtype=float)
+                # if type(surface) is not str:
+                #     for s in surface:
+                #         self.df.loc[idx, 'surface'] = s
+                # else:
+                #     self.df.loc[idx, 'surface'] = surface
                 self.df.loc[idx, 'surface'] = surface
 
                 self.calculate_gradient(idx)

gempy/core/data_modules/stack.py

@@ -215,7 +215,8 @@ def reset_order_series(self):
         """
         Reset the column order series to monotonic ascendant values.
         """
-        self.df.at[:, 'order_series'] = pn.RangeIndex(1, self.df.shape[0] + 1)
+        # self.df['order_series'] = pn.RangeIndex(1, self.df.shape[0] + 1) pandas 1.4 change
+        self.df.at[:, 'order_series'] = pn.RangeIndex(1, self.df.shape[0] + 1) # pandas 1.3 version
 
     @_setdoc_pro(reset_order_series.__doc__)
     def set_series_index(self, series_order: Union[list, np.ndarray], reset_order_series=True):

gempy/plot/visualization_2d.py

@@ -455,8 +455,8 @@ def plot_data(self, ax, section_name=None, cell_number=None, direction='y',
         points_df = points[select_projected_p]
         points_df['colors'] = points_df['surface'].map(self._color_lot)
 
-        points_df.plot.scatter(x=x, y=y, ax=ax, c='colors', s=70,  zorder=102,
-                               edgecolors='white',
+        points_df.plot.scatter(x=x, y=y, ax=ax, c=points_df['surface'].map(self._color_lot),
+                               s=70,  zorder=102, edgecolors='white',
                                colorbar=False)
         # points_df.plot.scatter(x=x, y=y, ax=ax, c='white', s=80,  zorder=101,
         #                        colorbar=False)

gempy/utils/docstring.py

@@ -10,7 +10,7 @@
 pole_vector = '(numpy.ndarray[float, 3]): 2D numpy array where axis 1 is the gradient values G_x, G_y, G_z of the pole while axis 0 is n number of' \
               ' orientations.'
 
-orientations = '(numpy.ndarray[float, 3]): 2D numpy array where axis 1 is are orientation values [dip, azimuth, polarity] of the pole while axis 0' \
+orientations = '(numpy.ndarray[float, 3]): 2D numpy array where axis 1 is are orientation values [azimuth, dip, polarity] of the pole while axis 0' \
                ' is n number of orientations. --- ' \
                '*Dip* is the inclination angle of 0 to 90 degrees measured from the horizontal plane downwards. ' \
                '*Azimuth* is the dip direction defined by a 360 degrees clockwise rotation, i.e. 0 = North,' \

gempy/utils/export.py

@@ -3,6 +3,7 @@
 from gempy.plot._visualization_2d import PlotData2D
 import numpy as np
 import os
+import itertools as it
 
 
 def export_geomap2geotiff(path, geo_model, geo_map=None, geotiff_filepath=None):
@@ -135,6 +136,147 @@ def export_moose_input(geo_model, path=None, filename='geo_model_units_moose_inp
 
     print("Successfully exported geological model as moose input to "+path)
 
+def export_shemat_suite_input_file(geo_model, path: str=None, filename: str='geo_model_SHEMAT_input'):
+    """
+    Method to export a 3D geological model as SHEMAT-Suite input-file for a conductive HT-simulation. 
+
+    Args:
+        path (str): Filepath for the exported input file (default './')
+        filename (str): name of exported input file (default 'geo_model_SHEMAT_input')
+    """
+    # get model dimensions
+    nx, ny, nz = geo_model.grid.regular_grid.resolution
+    xmin, xmax, ymin, ymax, zmin, zmax = geo_model.solutions.grid.regular_grid.extent
+    
+    delx = (xmax - xmin)/nx
+    dely = (ymax - ymin)/ny
+    delz = (zmax - zmin)/nz
+    
+    # get unit IDs and restructure them
+    ids = np.round(geo_model.solutions.lith_block)
+    ids = ids.astype(int)
+    
+    liths = ids.reshape((nx, ny, nz))
+    liths = liths.flatten('F')
+
+    # group litho in space-saving way
+    sequence = [len(list(group)) for key, group in it.groupby(liths)]
+    unit_id = [key for key, group in it.groupby(liths)]
+    combined = ["%s*%s" % (pair) for pair in zip(sequence,unit_id)]
+
+    combined_string = " ".join(combined)
+
+    # get number of units and set units string
+    units = geo_model.surfaces.df[['surface', 'id']]
+    
+    unitstring = ""
+    for index, rows in units.iterrows():
+        unitstring += f"0.01d-10    1.d0  1.d0  1.e-14	 1.e-10  1.d0  1.d0  3.74	0.  2077074.  10  2e-3	!{rows['surface']} \n" 	
+
+    # input file as f-string
+    fstring = f"""!==========>>>>> INFO
+# Title
+{filename}
+
+# linfo
+1 2 1 1
+
+# runmode
+1
+
+# timestep control
+0
+1           1           0           0
+
+# tunit
+1
+ 
+# time periods, records=1
+0      60000000    200      lin
+           
+# output times, records=10
+1
+6000000
+12000000
+18000000
+24000000
+30000000
+36000000
+42000000
+48000000
+54000000
+    
+# file output: hdf vtk
+
+# active temp
+
+# PROPS=bas
+
+# USER=none
+
+
+# grid
+{nx} {ny} {nz}
+
+# delx
+{nx}*{delx}
+
+# dely
+{ny}*{dely}
+
+# delz
+{nz}*{delz}
+
+!==========>>>>> NONLINEAR SOLVER
+# nlsolve
+50 0
+
+!==========>>>>> FLOW
+# lsolvef (linear solver control)
+1.d-12 64 500
+# nliterf (nonlinear iteration control)
+1.0d-10 1.0
+
+!==========>>>>> TEMPERATURE
+# lsolvet (linear solver control)
+1.d-12 64 500
+# nlitert (nonlinear iteration control)
+1.0d-10 1.0
+
+!==========>>>>> INITIAL VALUES
+# temp init
+{nx*ny*nz}*15.0d0  
+
+# head init
+{nx*ny*nz}*7500
+
+!==========>>>>> UNIT DESCRIPTION
+!!
+# units
+{unitstring}
+
+!==========>>>>>   define boundary properties
+# temp bcd, simple=top, value=init
+
+# temp bcn, simple=base, error=ignore
+{nx*ny}*0.06
+
+# uindex
+{combined_string}"""
+
+    if not path:
+        path = './'
+    if not os.path.exists(path):
+        os.makedirs(path)
+
+    f = open(path+filename, 'w+')
+    
+    f.write(fstring)
+    f.close()
+    
+    print("Successfully exported geological model as SHEMAT-Suite input to "+path) 
+
+
 def export_pflotran_input(geo_model, path=None, filename='pflotran.ugi'):
     """
     Method to export a 3D geological model as PFLOTRAN implicit unstructured grid
@@ -245,10 +387,11 @@ def export_flac3D_input(geo_model, path=None, filename='geomodel.f3grid'):
     vertices, elements, groups = __build_vertices_elements_groups__(geo_model)
 
     #open output file
-    if not path:
-        path = './'
-    if not os.path.exists(path):
-        os.makedirs(path)
+    #if not path:
+    #    path = './'
+    #if not os.path.exists(path):
+    #    os.makedirs(path)
+
     out = open(path+filename, 'w')
 
     #write gridpoints

setup.py

@@ -10,6 +10,7 @@
     packages=find_packages(exclude=('test', 'docs', 'examples')),
     include_package_data=True,
     install_requires=[
+        'pandas==1.3.4',
         'Theano>=1.0.4',
         'matplotlib',
         'numpy',