General bug fixes. Added option to clean only fractures touching the …

…boundary. and folder restructuring. Added the paper_materials folder

fracability/AbstractClasses.py

@@ -298,17 +298,17 @@ def save_shp(self, path: str):
                 for f_set in set_n:
                     final_path = os.path.join(output_path, f'{self.name}_{f_set}.shp')
                     entity_df = self.entity_df.loc[self.entity_df['f_set'] == f_set, :]
-                    entity_df.to_file(final_path, crs=self.crs)
+                    entity_df.to_file(final_path, crs=self.crs, engine='fiona')
             elif self.name == 'Boundary':
                 group_n = list(set(self.entity_df['b_group']))
                 for b_group in group_n:
                     final_path = os.path.join(output_path, f'{self.name}_{b_group}.shp')
                     entity_df = self.entity_df.loc[self.entity_df['b_group'] == b_group, :]
-                    entity_df.to_file(final_path, crs=self.crs)
+                    entity_df.to_file(final_path, crs=self.crs, engine='fiona')
             elif self.name == 'Nodes':
                 final_path = os.path.join(output_path, f'{self.name}.shp')
                 entity_df = self.entity_df
-                entity_df.to_file(final_path, crs=self.crs)
+                entity_df.to_file(final_path, crs=self.crs, engine='fiona')
 
             qgis_style_paths = QgisStyle().available_paths
 

fracability/Entities.py

@@ -1250,28 +1250,36 @@ def check_network(self, check_single=True, save_shp=None):
         else:
             print(overlaps_list_dict)
 
-    def clean_network(self, buffer = 0.05, inplace=True):
+    def clean_network(self, buffer = 0.05, inplace=True, only_boundary=False):
         """Tidy the intersection of the active entities in the fracture network. A buffer is applied to all the
         geometries to ensure intersection in a given radius.
 
+        :param only_boundary: Apply cleaning only on the fractures intersecting the boundary
         :param buffer: Applied buffer to the geometries of the entity.
         :param inplace: If true automatically replace the network with the clean one, if false then return the clean
          geopandas dataframe. Default is True
          """
 
         if inplace:
-            Geometry.tidy_intersections(self)
+            if only_boundary:
+                Geometry.tidy_intersections_boundary_only(self)
+            else:
+                Geometry.tidy_intersections(self)
         else:
-            return Geometry.tidy_intersections(self, inplace=False)
+            if only_boundary:
+                Geometry.tidy_intersections_boundary_only(self, inplace=False)
+            else:
+                Geometry.tidy_intersections(self, inplace=False)
 
-    def calculate_topology(self, clean_network=True):
+    def calculate_topology(self, clean_network=True, only_boundary=False):
         """
         Calculate the topology of the network and add the calculated nodes to the network.
 
         :param clean_network: If true, before calculating the topology the network is cleaned with the clean_network. Default is True
+        :param only_boundary: Apply cleaning only on the fractures intersecting the boundary
         """
         if clean_network is True:
-            self.clean_network()
+            self.clean_network(only_boundary=only_boundary)
 
         nodes_dict, origin_dict = Topology.nodes_conn(self)
         self.add_nodes_from_dict(nodes_dict,origin_dict=origin_dict, classes=None)

fracability/Plotters.py

@@ -73,7 +73,6 @@ def matplot_nodes(entity,
 
     1. I nodes: blue red circle
     2. Y nodes: green triangle
-    3. Y2 nodes: cyan triangle
     4. X nodes: blue square
     5. U nodes: yellow pentagon
     """
@@ -91,12 +90,10 @@ def matplot_nodes(entity,
     Y = np.where(node_types == 3)
     X = np.where(node_types == 4)
     U = np.where(node_types == 5)
-    Y2 = np.where(node_types == 6)
 
     # todo change the colors to reflect the pallete of the paper
     ax.plot(points[I][:, 0], points[I][:, 1], 'or', markersize=markersize)
     ax.plot(points[Y][:, 0], points[Y][:, 1], '^g', markersize=markersize)
-    ax.plot(points[Y2][:, 0], points[Y2][:, 1], '^c', markersize=markersize)
     ax.plot(points[X][:, 0], points[X][:, 1], 'sb', markersize=markersize)
     ax.plot(points[U][:, 0], points[U][:, 1], 'py', markersize=markersize)
 
@@ -381,12 +378,10 @@ def vtkplot_nodes(entity,
         3: 'Y',
         4: 'X',
         5: 'U',
-        6: 'Y2'
     }
     cmap_dict = {
         'I': 'Blue',
         'Y': 'Green',
-        'Y2': 'Cyan',
         'X': 'Red',
         'U': 'Yellow'
     }
@@ -407,7 +402,7 @@ def vtkplot_nodes(entity,
 
     actor = plotter.add_mesh(nodes,
                              scalars=class_names,
-                             render_points_as_spheres=True,
+                             render_points_as_spheres=False,
                              point_size=markersize,
                              show_scalar_bar=True,
                              scalar_bar_args=sargs,
@@ -559,10 +554,6 @@ def vtkplot_frac_net(entity,
     fractures = entity.fractures
     boundaries = entity.boundaries
 
-    if nodes is not None:
-        node_actor = vtkplot_nodes(nodes, markersize=markersize, return_plot=True)
-        plotter.add_actor(node_actor)
-
     if fractures is not None:
         fractures_actor = vtkplot_fractures(fractures, linewidth=fracture_linewidth,
                                             color=fracture_color, color_set=color_set, return_plot=True)
@@ -573,6 +564,11 @@ def vtkplot_frac_net(entity,
                                             color=boundary_color, return_plot=True)
         plotter.add_actor(boundary_actor)
 
+    if nodes is not None:
+        node_actor = vtkplot_nodes(nodes, markersize=markersize, return_plot=True)
+
+        plotter.add_actor(node_actor)
+
     if return_plot:
         actors = plotter.actors
         return actors

fracability/operations/Geometry.py

@@ -62,6 +62,51 @@ def tidy_intersections(obj, buffer=0.05, inplace: bool = True):
             line1 = gdf.loc[
                 idx_line1, 'geometry']  # Use as the reference line (in the int_node function) the new geometry.
 
+    print('\n\n')
+
+    if inplace:
+        obj.entity_df = gdf
+    else:
+        copy_obj = deepcopy(obj)
+        copy_obj.entity_df = gdf
+        return copy_obj
+
+def tidy_intersections_boundary_only(obj, buffer=0.05, inplace: bool = True):
+    """Method used to tidy shapefile intersections with the boundary of a fracture or fracture network object."""
+    if obj.name == 'FractureNetwork':
+        gdf = obj.fracture_network_to_components_df()
+        gdf = gdf.loc[gdf['type'] != 'node']
+    elif obj.name == 'Fractures':
+        gdf = obj.entity_df.copy()
+    else:
+        print('Cannot tidy intersection for nodes or only boundaries')
+        return
+    df_buffer = gdf.buffer(buffer)
+    print('\n\n')
+    for idx_line1, line in gdf.iterrows():
+        print(f'Calculating intersections on fracture: {idx_line1+1}/{len(gdf.index)}', end='\r')
+
+        if line['type'] == 'boundary':
+            continue
+
+        line1 = line['geometry']
+
+        idx_list = df_buffer.index[df_buffer.intersects(line1) == True]  # Subset the intersecting lines
+        idx_list = idx_list[idx_list != idx_line1]  # Exclude the reference line index to avoid self intersection
+
+        intersections = gdf.loc[idx_list]
+        for line2, idx_line2 in zip(intersections['geometry'], idx_list):  # for each intersecting line:
+            if intersections['type'][idx_line2] == 'boundary':
+                new_geom = int_node(line1, line2, [idx_line1, idx_line2], gdf)  # Calculate and add the intersection node.
+
+                for key, value in new_geom.items():
+                    gdf.loc[key, 'geometry'] = value  # substitute the original geometry with the new geometry
+
+                line1 = gdf.loc[
+                    idx_line1, 'geometry']  # Use as the reference line (in the int_node function) the new geometry.
+            else:
+                continue
+    print('\n\n')
     if inplace:
         obj.entity_df = gdf
     else:

paper_materials/Pontrelli/Pontrelli_analysis.py

@@ -0,0 +1,45 @@
+'''
+Script used to analyse the Pontrelli dataset.
+'''
+
+from fracability import Entities, Statistics  # import the Entities class
+
+data_path = 'Pontrelli_in/Set_a.shp'
+boundary_path = 'Pontrelli_in/Interpretation_boundary.shp'
+
+# Create the fractures and boundary objects.
+set_a = Entities.Fractures(shp=data_path, set_n=1)  # to add your data put the absolute path of the shp file
+
+boundary = Entities.Boundary(shp=boundary_path, group_n=1)
+
+fracture_net = Entities.FractureNetwork()
+
+fracture_net.add_fractures(set_a)
+
+fracture_net.add_boundaries(boundary)
+
+fracture_net.calculate_topology()
+
+fracture_net.vtk_plot(notebook=False)
+
+fracture_net.fractures.save_csv('Pontrelli_out')
+
+
+fitter = Statistics.NetworkFitter(fracture_net)
+
+fitter.fit('lognorm')
+fitter.fit('expon')
+fitter.fit('norm')
+fitter.fit('gamma')
+fitter.fit('powerlaw')
+fitter.fit('weibull_min')
+
+
+fitter.plot_PIT(bw=True)
+fitter.plot_summary(position=[1], sort_by='Mean_rank')
+# fitter.fit_result_to_clipboard()
+# Plotters.matplot_stats_ranks(fitter)
+
+
+
+

paper_materials/Pontrelli/Pontrelli_in/Interpretation_boundary.cpg

@@ -0,0 +1 @@
+UTF-8

paper_materials/Pontrelli/Pontrelli_in/Interpretation_boundary.prj

@@ -0,0 +1 @@
+PROJCS["WGS_1984_UTM_Zone_33N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",15.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]

paper_materials/Pontrelli/Pontrelli_in/Set_a.cpg

@@ -0,0 +1 @@
+UTF-8

paper_materials/Pontrelli/Pontrelli_in/Set_a.prj

@@ -0,0 +1 @@
+PROJCS["WGS_1984_UTM_Zone_33N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",15.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]