Merge branch 'master' of github.com:Loop3D/LoopStructural

Author: lachlangrose

LoopStructural/export/exporters.py

@@ -0,0 +1,175 @@
+"""
+Routines to export geological model data to file in a variety of formats
+"""
+import logging
+from pyevtk.hl import unstructuredGridToVTK, pointsToVTK
+from pyevtk.vtk import VtkTriangle
+import numpy as np
+
+from LoopStructural.utils.helper import create_box
+from LoopStructural.export.file_formats import FileFormat
+
+ 
+logger = logging.getLogger(__name__)
+
+
+def write_cubeface(model, file_name, data_label, nsteps, file_format):
+    """
+    Writes out the model as a cuboid with six rectangular surfaces
+
+    Parameters
+    ----------
+    model : GeologicalModel object
+        Geological model to export
+    file_name : string
+        Name of file that model is exported to, including path, but without the file extension
+    data_label : string
+        A data label to insert into export file
+    nsteps : np.array([num-x-steps, num-y-steps, num-z-steps])
+        3d array dimensions
+    file_format: export.fileformats.FileFormat object
+        Desired format of exported file
+
+    Returns
+    -------
+    True if successful
+
+    """
+    if file_format == FileFormat.VTK:
+        return _write_cubeface_evtk(model, file_name, data_label, nsteps)
+
+    logger.warning("Cannot export to file - format {} not supported yet".format(str(file_format)))
+    return False
+
+
+def write_vol(model, file_name, data_label, nsteps, file_format):
+    """
+    Writes out the model as a 3d volume grid
+
+    Parameters
+    ----------
+    model : GeologicalModel object
+        Geological model to export
+    file_name : string
+        Name of file that model is exported to, including path, but without the file extension
+    data_label : string
+        A data label to insert into export file
+    nsteps : np.array([num-x-steps, num-y-steps, num-z-steps])
+        3d array dimensions
+    file_format: export.fileformats.FileFormat object
+        Desired format of exported file
+
+    Returns
+    -------
+    True if successful
+
+    """
+    if file_format == FileFormat.VTK:
+        return _write_vol_evtk(model, file_name, data_label, nsteps)
+
+    logger.warning("Cannot export to file - format {} not supported yet".format(str(file_format)))
+    return False
+
+
+def _write_cubeface_evtk(model, file_name, data_label, nsteps, real_coords=True):
+    """
+    Writes out the model as a cuboid with six rectangular surfaces
+
+    Parameters
+    ----------
+    model : GeologicalModel object
+        Geological model to export
+    file_name : string
+        Name of file that model is exported to, including path, but without the file extension
+    data_label : string
+        A data label to insert into export file
+    nsteps : np.array([num-x-steps, num-y-steps, num-z-steps])
+        3d array dimensions
+
+    Returns
+    -------
+    True if successful
+
+    """
+    # Evaluate model at points
+    points, tri = create_box(model.bounding_box, nsteps)
+    val = model.evaluate_model(points, scale=False)
+    if real_coords:
+        model.rescale(points)
+
+    # Define vertices
+    x = np.zeros(points.shape[0])
+    y = np.zeros(points.shape[0])
+    z = np.zeros(points.shape[0])
+    for i in range(points.shape[0]):
+        x[i], y[i], z[i] = points[i][0], points[i][1], points[i][2]
+
+    # Define connectivity or vertices that belongs to each element
+    conn = np.zeros(tri.shape[0] * 3)
+    for i in range(tri.shape[0]):
+        conn[i*3], conn[i*3+1], conn[i*3+2] = tri[i][0], tri[i][1], tri[i][2]
+
+    # Define offset of last vertex of each element
+    offset = np.zeros(tri.shape[0])
+    for i in range(tri.shape[0]):
+        offset[i] = (i+1)*3
+
+    # Define cell types
+    ctype = np.full(tri.shape[0], VtkTriangle.tid)
+
+    try:
+        unstructuredGridToVTK(file_name, x, y, z, connectivity = conn, offsets = offset, cell_types = ctype, cellData = None, pointData = {data_label: val})
+    except Exception as e:
+        logger.warning("Cannot export cuboid surface to file {}: {}".format(file_name, str(e)))
+        return False
+    return True 
+
+
+def _write_vol_evtk(model, file_name, data_label, nsteps, real_coords=True):
+    """
+    Writes out the model as a 3d volume grid
+
+    Parameters
+    ----------
+    model : GeologicalModel object
+        Geological model to export
+    file_name : string
+        Name of file that model is exported to, including path, but without the file extension
+    data_label : string
+        A data label to insert into export file
+    nsteps : np.array([num-x-steps, num-y-steps, num-z-steps])
+        3d array dimensions
+
+    Returns
+    -------
+    True if successful
+
+    """
+    # Define grid spacing
+    loop_X = np.linspace(model.bounding_box[0, 0], model.bounding_box[1, 0], nsteps[0])
+    loop_Y = np.linspace(model.bounding_box[0, 1], model.bounding_box[1, 1], nsteps[1])
+    loop_Z = np.linspace(model.bounding_box[0, 2], model.bounding_box[1, 2], nsteps[2])
+
+    # Generate model values in 3d grid
+    xx, yy, zz = np.meshgrid(loop_X, loop_Y, loop_Z, indexing='ij')
+    xyz = np.array([xx.flatten(), yy.flatten(), zz.flatten()]).T
+    vals = model.evaluate_model(xyz, scale=False)
+    if real_coords:
+        model.rescale(xyz)
+
+    # Define vertices
+    x = np.zeros(xyz.shape[0])
+    y = np.zeros(xyz.shape[0])
+    z = np.zeros(xyz.shape[0])
+    for i in range(xyz.shape[0]):
+        x[i], y[i], z[i] = xyz[i][0], xyz[i][1], xyz[i][2]
+
+    # Write to grid
+    try:
+        pointsToVTK(file_name, x, y, z, data= { data_label: vals})
+    except Exception as e:
+        logger.warning("Cannot export volume to file {}: {}".format(file_name, str(e)))
+        return False
+    return True 
+
+

LoopStructural/export/file_formats.py

@@ -0,0 +1,12 @@
+"""
+Exported file formats
+"""
+from enum import Enum
+
+class FileFormat(Enum):
+    """ Enumeration of file export formats 
+    """
+    OBJ = 1
+    VTK = 2
+    GZ = 3 # Not supported yet
+    GLTF = 4 # Not supported yet

LoopStructural/interpolators/surfe_wrapper.py

@@ -76,7 +76,6 @@ def _setup_interpolator(self, **kwargs):
         if regression > 0:
             logger.info("Using regression smoothing %f" % regression)
             self.surfe.SetRegressionSmoothing(True,regression)
-            print("adding regression {}".format(regression))
         greedy = kwargs.get("greedy", (0, 0))
 
         if greedy[0] > 0 or greedy[1] > 0:

LoopStructural/modelling/features/structural_frame.py

@@ -14,8 +14,7 @@ class StructuralFrame:
     def __init__(self, name, features, fold=None):
         """
         Structural frame is a curvilinear coordinate system defined by
-        structural
-        observations associated with a fault or fold.
+        structural observations associated with a fault or fold.
 
         Parameters
         ----------
@@ -52,7 +51,7 @@ def set_model(self, model):
         self.model = model
         for f in self.features:
             if f is None:
-                continue
+                continue    
             f.set_model(model)
 
     def add_region(self, region):

LoopStructural/visualisation/model_visualisation.py

@@ -13,7 +13,7 @@
     logger.error("Please install lavavu: pip install lavavu")
 import numpy as np
 from skimage.measure import marching_cubes_lewiner as marching_cubes
-
+from LoopStructural.modelling.features import GeologicalFeature
 from LoopStructural.utils.helper import create_surface, get_vectors, create_box
 
 # adapted/copied from pyvista for sphinx scraper
@@ -146,7 +146,7 @@ def add_section(self, geological_feature=None, axis='x', value=None, **kwargs):
         -------
 
         """
-
+        print('aa')
         if axis == 'x':
             tri, yy, zz = create_surface(self.bounding_box[:, [1, 2]], self.nsteps[[1, 2]])
             xx = np.zeros(zz.shape)
@@ -178,12 +178,13 @@ def add_section(self, geological_feature=None, axis='x', value=None, **kwargs):
         points[:, 2] = zz
 
         surf = self.lv.triangles(name)
-        surf.vertices(self.model.rescale(points))
+        surf.vertices(self.model.rescale(points,inplace=False))
         surf.indices(tri)
         logger.info("Adding %s section at %f" % (axis, value))
         if geological_feature is None:
             surf.colours(colour)
-        if geological_feature is not None:
+
+        if geological_feature is not None and type(geological_feature) == GeologicalFeature:
             if 'norm' in kwargs:
                 surf.values(np.linalg.norm(
                     geological_feature.evaluate_gradient(points), axis=1),
@@ -198,7 +199,17 @@ def add_section(self, geological_feature=None, axis='x', value=None, **kwargs):
             logger.info("Colouring section with %s min: %f, max: %f" % (
                 geological_feature.name, geological_feature.min(), geological_feature.max()))
             surf.colourmap(cmap, range=[geological_feature.min(), geological_feature.max()])
-
+        if geological_feature == 'model' and self.model is not None:
+            name = kwargs.get('name','model_section')
+            surf.values(self.model.evaluate_model(points,scale=True),
+                            name)
+            surf["colourby"] = name
+            cmap = lavavu.cubehelix(100)
+            if 'cmap' in kwargs:
+                cmap = kwargs['cmap']
+            # logger.info("Colouring section with %s min: %f, max: %f" % (
+            #     geological_feature.name, geological_feature.min(), geological_feature.max()))
+            # surf.colourmap(cmap, range=[geological_feature.min(), geological_feature.max()])
     def add_isosurface(self, geological_feature, value = None, isovalue=None,
                      paint_with=None, slices=None, colour='red', nslices=None, 
                      cmap=None, filename=None, names=None, colours=None,**kwargs):

README.md

@@ -16,7 +16,7 @@ Loop is led by Laurent Ailleres (Monash University) with a team of Work Package
 
 The project benefits from in-kind contributions from the Geological Survey of Canada, the British Geological Survey, the French Bureau de Recherches Geologiques et Minieres, the RING group at the Universite de Lorraine, the RWTH in Aachen, Germany and AUSCOPE
 
-* Python/cython implementation of a piecewise linear interpolation (Discrete Smooth Interpolator in Gocad) 
+* Python/cython implementation of a Discrete interpolatiors  
 * Fold interpolation using constraints outlined in Laurent 2016 with fold geostatistical tools shown in Grose et al., 2017
 * Fault interpolation 
 

conda/meta.yaml

@@ -1,5 +1,5 @@
 {% set name = "loopstructural" %}
-{% set version = "1.0.4" %}
+{% set version = "1.0.5" %}
 
 package:
   name: "{{ name|lower }}"
@@ -14,20 +14,22 @@ build:
   script: "{{ PYTHON }} -m pip install ."
 
 requirements:
+  build:
+    - cython
   host:
     - pip
     - python
     - cython
     - setuptools
     - numpy
-  run:
-    - Cython
-    - numpy
-    - pandas
-    - scipy
-    - matplotlib
-    - scikit-image
-#     - scikit-learn
+  # run:
+  #   - Cython
+  #   - numpy
+  #   - pandas
+  #   - scipy
+  #   - matplotlib
+  #   - scikit-image
+  #   - scikit-learn
 # test:
 #   imports:
 #     - LoopStructural

requirements.txt

@@ -6,5 +6,6 @@ numpy
 matplotlib
 setuptools
 pandas
+pyevtk
 
 

setup.py

@@ -20,14 +20,14 @@ def get_version(rel_path):
 setup(
 	name="LoopStructural",
 	install_requires=[
-	'Cython',
-	'numpy',
-	'pandas',
-	'scipy',
-	'matplotlib',
-	# 'lavavu',
-	'scikit-image',
-	'scikit-learn'
+	# 'Cython',
+	# 'numpy',
+	# 'pandas',
+	# 'scipy',
+	# 'matplotlib',
+	# # 'lavavu',
+	# 'scikit-image',
+	# 'scikit-learn'
 	],
 	version=get_version("LoopStructural/__init__.py"),
     packages=find_packages(),

tests/unit_tests/modelling/test_fold_event.py

@@ -1 +1,11 @@
-from LoopStructural.modelling.fold import FoldEvent
+from LoopStructural.modelling.fold import FoldEvent
+from LoopStructural.modelling.features import StructuralFrame
+def test_constructor():
+    
+    pass
+
+def test_constant_fold_axis():
+    pass
+
+def test_rotation_fold_axis():
+    pass