Skip to content

Commit

Permalink
changes for Adapters.py. Finalized the spacing script example
Browse files Browse the repository at this point in the history
  • Loading branch information
gbene committed Nov 24, 2023
1 parent 67bf65f commit 7e17a78
Show file tree
Hide file tree
Showing 5 changed files with 39 additions and 105 deletions.
Binary file modified development/Spacing/data/cava_pontrelli/FN_set_1.dbf
Binary file not shown.
Binary file modified development/Spacing/data/cava_pontrelli/FN_set_1.shp
Binary file not shown.
Binary file modified development/Spacing/data/cava_pontrelli/FN_set_1.shx
Binary file not shown.
112 changes: 26 additions & 86 deletions development/Spacing/main.py
Original file line number Diff line number Diff line change
Expand Up @@ -47,31 +47,6 @@ def centers_to_lines(center_coords, lengths, frac_dir, assign_id=True) -> pv.Pol
return lines


def plot_rose(dir_data):
bin_edges = np.arange(-5, 366, 10)
number_of_strikes, bin_edges = np.histogram(dir_data, bin_edges)
number_of_strikes[0] += number_of_strikes[-1]
half = np.sum(np.split(number_of_strikes[:-1], 2), 0)
two_halves = np.concatenate([half, half])
fig = plt.figure(figsize=(16, 8))

ax = fig.add_subplot(111, projection='polar')

ax.bar(np.deg2rad(np.arange(0, 360, 10)), two_halves,
width=np.deg2rad(10), bottom=0.0, color='.8', edgecolor='k')
ax.set_theta_zero_location('N')
ax.set_theta_direction(-1)
ax.set_title('Rose Diagram of the "Fault System"', y=1.10, fontsize=15)

fig.tight_layout()
plt.show()


def histogram(dir_data):
sns.histplot(dir_data, bins=18)
plt.show()


def center_points(center_coords, lengths, dir, resolution=10):
"""Function used to create a series of points along a line with center, length and direction."""

Expand Down Expand Up @@ -99,9 +74,9 @@ def center_points(center_coords, lengths, dir, resolution=10):
df = gpd.read_file(data_path)
boundary_df = gpd.read_file(boundary_path)

length_plane = 120
length_scanl = 150
resolution = 10
length_plane = 150
length_scanl = 170
resolution = 100


df_faults = df.loc[df['Fault'] == 1]
Expand Down Expand Up @@ -137,72 +112,37 @@ def center_points(center_coords, lengths, dir, resolution=10):
fractures_diss = fractures.entity_df.dissolve()

split_scanlines_list = []

# for scanline in scanlines_shp:
# split_geom = split(scanline, boundary_shp)
# for segment in split_geom.geoms:
# if boundary_shp.buffer(0.001).contains(segment):
# split_scanlines_list.append(segment)

for scanline in scanlines_shp:
split_geom = split(scanline, boundary_shp)
split_geom = split(scanline, fractures_diss['geometry'].values[0])
for segment in split_geom.geoms:
if boundary_shp.covers(segment):
split_segment = split(segment, fractures_diss['geometry'].values[0])
split_scanlines_list.append(list(split_segment.geoms))
split_segment = split(segment, boundary_shp)
for line in split_segment.geoms:
if boundary_shp.buffer(0.001).contains(line):
split_scanlines_list.append(line)


split_scanlines_df = gpd.GeoDataFrame(geometry=np.concatenate(split_scanlines_list))
split_scanlines_df = gpd.GeoDataFrame(geometry=split_scanlines_list)

valid_scanlines = Entities.Fractures(gdf=split_scanlines_df,set_n=1)
valid_scanlines = Entities.Fractures(gdf=split_scanlines_df, set_n=1)

frac_net = Entities.FractureNetwork()
# frac_net.add_fractures(valid_scanlines)
frac_net.add_fractures(scanlines)
frac_net.add_fractures(valid_scanlines)
# frac_net.add_fractures(scanlines)
frac_net.add_boundaries(boundary)

frac_net.fractures.vtk_plot()

# frac_net.calculate_topology()
#
#
# fitter = Statistics.NetworkFitter(frac_net)
# fitter.fit('lognorm')
# matplot_stats_summary(fitter.get_fitted_distribution('lognorm'))

#
#
# fractures_diss = fractures.entity_df.dissolve()
# splitted_list = []
#
# for i, scanline in enumerate(scanlines.entity_df['geometry'].values):
#
# splitted = split(scanline, fractures_diss['geometry'].values[0])
# splitted_list.append(list(splitted.geoms))
#
# scanline_df = gpd.GeoDataFrame(geometry=np.concatenate(splitted_list))
#
# new_scanlines = Entities.Fractures(gdf=scanline_df, set_n=4)
#
#
#
# frac_net = Entities.FractureNetwork()
# frac_net.add_fractures(new_scanlines)
# frac_net.add_boundaries(boundary)
# print(frac_net.sets)
# frac_net.cut_net_on_boundary()
print(len(frac_net.fractures.entity_df['length']))
frac_net.calculate_topology()
# frac_net.vtk_plot()

# frac_net.calculate_topology()
# sns.histplot(new_scanlines.entity_df['length'].values)
# plt.show()


# fitter = Statistics.NetworkFitter(frac_net)
# fitter.fit('lognorm')
# matplot_stats_summary(fitter.get_fitted_distribution('lognorm'))

# print(new_scanlines.vtk_object.array_names)
# print(frac_net.fraction_censored)
#
# plotter = pv.Plotter()
# plotter.background_color = 'gray'
# plotter.add_mesh(boundary.vtk_object, color='white')
# # plotter.add_mesh(line, color='r')
# # # plotter.add_mesh(mean_plane, color='b')
# # # plotter.add_points(test_points, color='y')
# plotter.add_mesh(lines)
# plotter.view_xy()
# plotter.enable_image_style()
# plotter.show()
fitter = Statistics.NetworkFitter(frac_net)
fitter.fit('lognorm')
matplot_stats_summary(fitter.get_fitted_distribution('lognorm'))
32 changes: 13 additions & 19 deletions src/fracability/Adapters.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,16 +5,15 @@
TODO:
+ Add a adapter abstract class in the AbstractClasses.
+ Use as input gpd for all
+ Use multiblock instead of appender for the fracture network
+ Try and find a way to crate vtk object directly from polygon and not from lines of polygon
+ Use multiblock instead of appender for fracture network (or all?)?
+ Try and find a way to crate boundary vtk object directly from polygon and not from lines of polygon
"""

import geopandas
import networkx
import numpy as np
from pyvista import PolyData, lines_from_points
from pyvista import PolyData, lines_from_points, MultiBlock

from vtkmodules.vtkFiltersCore import vtkAppendPolyData

Expand Down Expand Up @@ -50,21 +49,19 @@ def frac_vtk_rep(input_df: geopandas.GeoDataFrame) -> PolyData:

points = np.stack((x, y, z), axis=1) # Stack the coordinates to a [n,3] shaped array
# offset = np.round(points[0][0])
pv_obj = lines_from_points(points) # Create the corresponding vtk line with the given points
pv_obj.cell_data['type'] = ['fracture'] * pv_obj.GetNumberOfCells()
pv_obj.point_data['type'] = ['fracture'] * pv_obj.GetNumberOfPoints()

pv_obj.cell_data['f_set'] = [set_n] * pv_obj.GetNumberOfCells()
pv_obj.point_data['f_set'] = [set_n] * pv_obj.GetNumberOfPoints()
if len(points) != 0: # Avoid empty geometries
pv_obj = lines_from_points(points) # Create the corresponding vtk line with the given points
pv_obj.cell_data['type'] = ['fracture'] * pv_obj.GetNumberOfCells()
pv_obj.point_data['type'] = ['fracture'] * pv_obj.GetNumberOfPoints()

pv_obj.cell_data['RegionId'] = [index] * pv_obj.GetNumberOfCells()
pv_obj.point_data['RegionId'] = [index] * pv_obj.GetNumberOfPoints()
pv_obj.cell_data['f_set'] = [set_n] * pv_obj.GetNumberOfCells()
pv_obj.point_data['f_set'] = [set_n] * pv_obj.GetNumberOfPoints()

if 'length' in input_df.columns:
pv_obj.cell_data['length'] = input_df.loc[index, 'length']
pv_obj.cell_data['RegionId'] = [index] * pv_obj.GetNumberOfCells()
pv_obj.point_data['RegionId'] = [index] * pv_obj.GetNumberOfPoints()

# pv_obj.cell_data_to_point_data()
# line.plot()
if 'length' in input_df.columns:
pv_obj.cell_data['length'] = input_df.loc[index, 'length']

appender.AddInputData(pv_obj) # Add the new object

Expand Down Expand Up @@ -159,9 +156,6 @@ def networkx_rep(input_object: PolyData) -> networkx.Graph():
lines = np.delete(lines,
np.arange(0, lines.size, 3)).reshape(-1, 2) # remove padding eg. [2 id1 id2 2 id3 id4 ...] -> remove the 2

# test_types = np.array([{'type': t} for t in network['type']])
# edges = np.c_[lines.reshape(-1, 2)]

network = nx.Graph() # Create a networkx graph instance

network.add_edges_from(lines) # Add the edges
Expand Down

0 comments on commit 7e17a78

Please sign in to comment.