Merge pull request #29 from Loop3D/map2loop

Map2loop

Author: lachlangrose

.github/workflows/python-publish.yml

@@ -5,7 +5,7 @@ on:
     types: [created]
 
 jobs:
-  deploy:
+  flake8:
     runs-on: ubuntu-latest
     steps:
     - uses: actions/checkout@v2
@@ -16,12 +16,25 @@ jobs:
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
-        pip install twine flake8
+        pip install flake8
     - name: Lint with flake8 for syntax errors
       run: |
         pip install flake8
         flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
         flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
+  manylinux:
+    runs-on: ubuntu-latest
+    needs: flake8
+    steps:
+    - uses: actions/checkout@v2
+    - name: Set up Python
+      uses: actions/setup-python@v1
+      with:
+        python-version: 3.8
+    - name: Install dependencies
+      run: |
+        python -m pip install --upgrade pip
+        pip install twine 
     - name: Build manylinux Python wheels
       uses: RalfG/python-wheels-manylinux-build@v0.2.2-manylinux2010_x86_64
       with:
@@ -36,3 +49,37 @@ jobs:
         TWINE_PASSWORD: ${{ secrets.PYPI_PASSWORD }}
       run: |
         twine upload wheelhouse/*-manylinux*.whl
+        
+  build-windows:
+    runs-on: windows-latest
+    needs: flake8
+    strategy:
+      matrix:
+        python: ['3.6','3.7','3.8']
+    steps:
+    - uses: actions/checkout@v2
+    - uses: goanpeca/action-setup-conda@v1
+      with:
+        python-version: ${{ matrix.python }} 
+        activate-environment: loop    
+    - name: Installing dependencies
+      shell: bash -l {0}
+      run: |
+        python --version
+        pip install -r requirements.txt
+        conda info
+        conda list
+    - name: Building Loop wheel and installing
+      shell: bash -l {0}
+      run: |
+        python setup.py bdist_wheel
+        python setup.py bdist
+    - name: Publish wheels to PyPI
+      env:
+        TWINE_USERNAME: ${{ secrets.PYPI_USERNAME }}
+        TWINE_PASSWORD: ${{ secrets.PYPI_PASSWORD }}
+      shell: bash -l {0}
+      run : |
+        pip install twine
+        twine upload dist/*
+        

LoopStructural/modelling/core/geological_model.py

@@ -128,6 +128,31 @@ def from_map2loop_directory(cls, m2l_directory,**kwargs):
         m2l_data = process_map2loop(m2l_directory)
         return build_model(m2l_data,**kwargs), m2l_data
 
+    @classmethod
+    def from_file(cls, file):
+        try:
+            import dill as pickle
+        except ImportError:
+            logger.error("Cannot import from file, dill not installed")
+            return None
+        model = pickle.load(open(file,'rb'))
+        if type(model) == GeologicalModel:
+            return model
+        else:
+            logger.error('{} does not contain a geological model'.format(file))
+            return None
+
+    def to_file(self, file):
+        try:
+            import dill as pickle
+        except ImportError:
+            logger.error("Cannot write to file, dill not installed")
+            return
+        try:
+            pickle.dump(self,open(file,'wb'))
+        except pickle.PicklingError:
+            logger.error('Error saving file')
+
     def _add_feature(self, feature):
         """
         Add a feature to the model stack
@@ -262,7 +287,7 @@ def create_from_feature_list(self, features):
             if featuretype == 'folded_strati':
                 self.create_and_add_folded_foliation(f)
 
-    def get_interpolator(self, interpolatortype='PLI', nelements=5e5,
+    def get_interpolator(self, interpolatortype='PLI', nelements=1e5,
                          buffer=0.2, **kwargs):
         """
         Returns an interpolator given the arguments, also constructs a
@@ -622,6 +647,8 @@ def _add_domain_fault_above(self, feature):
 
         """
         for f in reversed(self.features):
+            if f.name == feature.name:
+                continue
             if f.type == 'domain_fault':
                 feature.add_region(lambda pos: f.evaluate_value(pos) < 0)
                 break
@@ -642,6 +669,8 @@ def _add_domain_fault_below(self, domain_fault):
 
         """
         for f in reversed(self.features):
+            if f.name == domain_fault.name:
+                continue
             f.add_region(lambda pos: domain_fault.evaluate_value(pos) > 0)
             if f.type == 'unconformity':
                 break
@@ -818,11 +847,14 @@ def create_and_add_domain_fault(self, fault_surface_data, **kwargs):
         # build feature
         domain_fault = domain_fault_feature_builder.build(**kwargs)
         domain_fault.type = 'domain_fault'
+        self._add_feature(domain_fault)
+        self._add_domain_fault_below(domain_fault)
+
         # uc_feature = UnconformityFeature(uc_feature_base,0)
         # iterate over existing features and add the unconformity as a
         # region so the feature is only
         # evaluated where the unconformity is positive
-        return self.add_unconformity(domain_fault, 0)
+        return domain_fault
 
     def create_and_add_fault(self, fault_surface_data, displacement, **kwargs):
         """
@@ -988,7 +1020,7 @@ def voxet(self, nsteps=(50, 50, 25)):
         """
         return {'bounding_box': self.bounding_box, 'nsteps': nsteps}
 
-    def regular_grid(self, nsteps=(50, 50, 25), shuffle = True):
+    def regular_grid(self, nsteps=(50, 50, 25), shuffle = True, rescale=True):
         """
         Return a regular grid within the model bounding box
 
@@ -1012,6 +1044,8 @@ def regular_grid(self, nsteps=(50, 50, 25), shuffle = True):
         locs = np.array([xx.flatten(), yy.flatten(), zz.flatten()]).T
         if shuffle:
             np.random.shuffle(locs)
+        if rescale:
+            locs = self.rescale(locs)
         return locs
 
     def evaluate_model(self, xyz, rescale=True):

LoopStructural/utils/map2loop.py

@@ -43,6 +43,8 @@ def process_map2loop(m2l_directory, flags={}):
     except:
         for g in groups['group'].unique():
             supergroups[g] = g
+    supergroups.pop('\n')
+
 
 
     bb = pd.read_csv(m2l_directory+'/tmp/bbox.csv')
@@ -85,6 +87,23 @@ def process_map2loop(m2l_directory, flags={}):
                 unit_id += 1
                 strat_val[c] = val[g]
                 val[g] += thickness[c]
+    group_name = None
+    for g, i in stratigraphic_column.items():
+        if len(i) ==0:
+            for gr, sg in supergroups.items():
+                if sg == g:
+                    group_name = gr
+                    break
+            try:
+                if group_name is None:
+                    continue
+                c=groups.loc[groups['group']==group_name,'code'].to_numpy()[0]
+                strat_val[c] = 0
+                stratigraphic_column[g] = {c:{'min':0,'max':9999,'id':unit_id}}
+                unit_id+=1
+                group_name = None
+            except:
+                print('Couldnt process {}'.format(g))
     contacts['val'] = np.nan
     for o in strat_val:
         contacts.loc[contacts['formation'] == o, 'val'] = strat_val[o]
@@ -140,7 +159,7 @@ def process_map2loop(m2l_directory, flags={}):
             'bounding_box':bb,
             'strat_va':strat_val}
 
-def build_model(m2l_data, skip_faults = False, fault_params = None, foliation_params=None):
+def build_model(m2l_data, skip_faults = False, unconformities=False, fault_params = None, foliation_params=None):
     """[summary]
 
     [extended_summary]
@@ -199,12 +218,12 @@ def build_model(m2l_data, skip_faults = False, fault_params = None, foliation_pa
 
     ## loop through all of the groups and add them to the model in youngest to oldest.
     group_features = []
-    for i in m2l_data['groups']['group number'].unique():
+    for i in np.sort(m2l_data['groups']['group number'].unique()):
         g = m2l_data['groups'].loc[m2l_data['groups']['group number'] == i, 'group'].unique()[0]
         group_features.append(model.create_and_add_foliation(g,
                                                             **foliation_params))
         # if the group was successfully added (not null) then lets add the base (0 to be unconformity)
-        # if group_features[-1]:
-        #     model.add_unconformity(group_features[-1], 0)
+        if group_features[-1] and unconformities:
+            model.add_unconformity(group_features[-1], 0)
     model.set_stratigraphic_column(m2l_data['stratigraphic_column'])
     return model

LoopStructural/visualisation/model_visualisation.py

@@ -149,8 +149,9 @@ def add_section(self, geological_feature=None, axis='x', value=None, **kwargs):
                 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, **kwargs):
+    def add_isosurface(self, geological_feature, value = None, isovalue=None,
+                     paint_with=None, slices=None, colour='red', nslices=None, 
+                     cmap=None, filename=None, **kwargs):
         """ Plot the surface of a geological feature 
 
         [extended_summary]
@@ -173,6 +174,8 @@ def add_isosurface(self, geological_feature, value = None, isovalue=None, paint_
             [description], by default None
         cmap : [type], optional
             [description], by default None
+        filename: string, optional
+            filename for exporting 
 
         Returns
         -------
@@ -188,7 +191,7 @@ def add_isosurface(self, geological_feature, value = None, isovalue=None, paint_
         # do isosurfacing of support using marching tetras/cubes
         x = np.linspace(self.bounding_box[0, 0], self.bounding_box[1, 0], self.nsteps[0])
         y = np.linspace(self.bounding_box[0, 1], self.bounding_box[1, 1], self.nsteps[1])
-        z = np.linspace(self.bounding_box[1, 2], self.bounding_box[0, 2], self.nsteps[2])
+        z = np.linspace(self.bounding_box[0, 2], self.bounding_box[1, 2], self.nsteps[2])
         xx, yy, zz = np.meshgrid(x, y, z, indexing='ij')
         points = np.array([xx.flatten(), yy.flatten(), zz.flatten()]).T
         val = geological_feature.evaluate_value(points)
@@ -240,9 +243,20 @@ def add_isosurface(self, geological_feature, value = None, isovalue=None, paint_
             except ValueError:
                 logger.warning("no surface to mesh, skipping")
                 continue
+            
+
             name = geological_feature.name
             name = kwargs.get('name', name)
             name += '_iso_%f' % isovalue
+            if filename is not None:
+                try:
+                    import meshio
+                except ImportError:
+                    logger.error("Could not save surfaces, meshio is not installed")
+                meshio.write_points_cells(filename.format(name),
+                self.model.rescale(verts),
+                [("triangle", faces)]
+                )
             surf = self.lv.triangles(name)
             surf.vertices(verts)
             surf.indices(faces)
@@ -353,12 +367,12 @@ def add_model_surfaces(self, faults = True, cmap='tab20', **kwargs):
             if g in self.model.feature_name_index:
                 feature = self.model.features[self.model.feature_name_index[g]]
                 for u, vals in self.model.stratigraphic_column[g].items():
-                    self.add_isosurface(feature, isovalue=vals['max'],name=u,colour=tab.colors[ci,:])
+                    self.add_isosurface(feature, isovalue=vals['max'],name=u,colour=tab.colors[ci,:],**kwargs)
                     ci+=1
         if faults:
             for f in self.model.features:
                 if f.type == 'fault':
-                    self.add_isosurface(f,isovalue=0)
+                    self.add_isosurface(f,isovalue=0,**kwargs)
 
 
     def add_vector_field(self, geological_feature, **kwargs):

examples/1_basic/plot_1_data_prepration.py

@@ -123,6 +123,14 @@
 # Geological feature can be evaluated:
 # * for the scalar field value at a location
 # * for the gradient of the scalar field at a location
+# To evaluate a model feature (scalar value or gradient) use the:
+# :code:`model.evaluate_feature_value(feature_name, locations)` or 
+# :code:`model.evaluate_feature_gradient(feature_name, locations)`
+# Where the feature_name is the string naming the feature and locations is a numpy array of 
+# xyz coordinates.
+# 
+# In the following example we will use matplotlib to visualise these results however, the
+# next tutorial will show how to use the lavavu visualisation model.  
 
 
 import matplotlib.pyplot as plt
@@ -135,10 +143,10 @@
 xx = np.zeros_like(yy)
 xx[:] = 5
 
-vals = conformable_feature.evaluate_value(model.scale(np.array([xx.flatten(),yy.flatten(),zz.flatten()]).T))
+vals = model.evaluate_feature_value('conformable',np.array([xx.flatten(),yy.flatten(),zz.flatten()]).T)
 fig, ax = plt.subplots(1,2,figsize=(20,10))
-ax[0].contourf(vals.reshape((100,100)))
-ax[0].contour(vals.reshape((100,100)),[0,1])
+ax[0].contourf(vals.reshape((100,100)),extent=(0,10,0,10))
+ax[0].contour(vals.reshape((100,100)),[0,1],extent=(0,10,0,10))
 
 # Y section
 x = np.linspace(0,10,100)
@@ -148,8 +156,8 @@
 yy = np.zeros_like(xx)
 yy[:] = 5
 
-vals = conformable_feature.evaluate_value(model.scale(np.array([xx.flatten(),yy.flatten(),zz.flatten()]).T))
-ax[1].contourf(vals.reshape((100,100)))
-ax[1].contour(vals.reshape((100,100)),[0,1])
+vals = model.evaluate_feature_value('conformable',np.array([xx.flatten(),yy.flatten(),zz.flatten()]).T)
+ax[1].contourf(vals.reshape((100,100)),extent=(0,10,0,10))
+ax[1].contour(vals.reshape((100,100)),[0,1],extent=(0,10,0,10))
 
 plt.show()