From 6f901ff9db22e3acee3d9ea9be7af929dacda959 Mon Sep 17 00:00:00 2001 From: rabii-chaarani Date: Mon, 28 Jul 2025 13:17:16 +0930 Subject: [PATCH 1/5] refactor: remove MapData dependency from thickness calculators --- map2loop/thickness_calculator.py | 68 +++++++++++++++++++++----------- 1 file changed, 44 insertions(+), 24 deletions(-) diff --git a/map2loop/thickness_calculator.py b/map2loop/thickness_calculator.py index 3da0ad40..f074f9b3 100644 --- a/map2loop/thickness_calculator.py +++ b/map2loop/thickness_calculator.py @@ -5,7 +5,8 @@ calculate_endpoints, multiline_to_line, find_segment_strike_from_pt, - set_z_values_from_raster_df + set_z_values_from_raster_df, + value_from_raster ) from .m2l_enums import Datatype from .interpolators import DipDipDirectionInterpolator @@ -16,6 +17,7 @@ # external imports from abc import ABC, abstractmethod +from typing import Optional import scipy.interpolate import beartype import numpy @@ -24,6 +26,7 @@ import geopandas import shapely import math +from osgeo import gdal from shapely.errors import UnsupportedGEOSVersionError class ThicknessCalculator(ABC): @@ -34,12 +37,19 @@ class ThicknessCalculator(ABC): ABC (ABC): Derived from Abstract Base Class """ - def __init__(self, max_line_length: float = None): + def __init__( + self, + dtm_data: Optional[gdal.Dataset] = None, + bounding_box: Optional[dict] = None, + max_line_length: Optional[float] = None + ): """ Initialiser of for ThicknessCalculator """ self.thickness_calculator_label = "ThicknessCalculatorBaseClass" self.max_line_length = max_line_length + self.dtm_data = dtm_data + self.bounding_box = bounding_box def type(self): """ @@ -58,7 +68,8 @@ def compute( stratigraphic_order: list, basal_contacts: geopandas.GeoDataFrame, structure_data: pandas.DataFrame, - map_data: MapData, + geology_data: geopandas.GeoDataFrame, + sampled_contacts: geopandas.GeoDataFrame, ) -> pandas.DataFrame: """ Execute thickness calculator method (abstract method) @@ -115,7 +126,8 @@ def compute( stratigraphic_order: list, basal_contacts: geopandas.GeoDataFrame, structure_data: pandas.DataFrame, - map_data: MapData, + geology_data: geopandas.GeoDataFrame, + sampled_contacts: geopandas.GeoDataFrame, ) -> pandas.DataFrame: """ Execute thickness calculator method takes unit data, basal_contacts and stratigraphic order and attempts to estimate unit thickness. @@ -144,7 +156,7 @@ def compute( basal_unit_list = basal_contacts["basal_unit"].to_list() sampled_basal_contacts = rebuild_sampled_basal_contacts( - basal_contacts=basal_contacts, sampled_contacts=map_data.sampled_contacts + basal_contacts=basal_contacts, sampled_contacts=sampled_contacts ) if len(stratigraphic_order) < 3: @@ -207,11 +219,16 @@ class InterpolatedStructure(ThicknessCalculator): -> pandas.DataFrame: Calculates a thickness map for the overall map area. """ - def __init__(self, max_line_length: float = None): + def __init__( + self, + dtm_data: Optional[gdal.Dataset] = None, + bounding_box: Optional[dict] = None, + max_line_length: Optional[float] = None + ): """ Initialiser for interpolated structure version of the thickness calculator """ - super().__init__(max_line_length) + super().__init__(dtm_data, bounding_box, max_line_length) self.thickness_calculator_label = "InterpolatedStructure" self.lines = None @@ -223,7 +240,8 @@ def compute( stratigraphic_order: list, basal_contacts: geopandas.GeoDataFrame, structure_data: pandas.DataFrame, - map_data: MapData, + geology_data: geopandas.GeoDataFrame, + sampled_contacts: geopandas.GeoDataFrame, ) -> pandas.DataFrame: """ Execute thickness calculator method takes unit data, basal_contacts, stratigraphic order, orientation data and @@ -264,7 +282,7 @@ def compute( # increase buffer around basal contacts to ensure that the points are included as intersections basal_contacts["geometry"] = basal_contacts["geometry"].buffer(0.01) # get the sampled contacts - contacts = geopandas.GeoDataFrame(map_data.sampled_contacts) + contacts = geopandas.GeoDataFrame(sampled_contacts) # build points from x and y coordinates geometry2 = geopandas.points_from_xy(contacts['X'], contacts['Y']) contacts.set_geometry(geometry2, inplace=True) @@ -272,8 +290,7 @@ def compute( # set the crs of the contacts to the crs of the units contacts = contacts.set_crs(crs=basal_contacts.crs) # get the elevation Z of the contacts - dtm_data = map_data.get_map_data(Datatype.DTM) - contacts = set_z_values_from_raster_df(dtm_data, contacts) + contacts = set_z_values_from_raster_df(self.dtm_data, contacts) # update the geometry of the contact points to include the Z value contacts["geometry"] = contacts.apply( lambda row: shapely.geometry.Point(row.geometry.x, row.geometry.y, row["Z"]), axis=1 @@ -281,12 +298,10 @@ def compute( # spatial join the contact points with the basal contacts to get the unit for each contact point contacts = contacts.sjoin(basal_contacts, how="inner", predicate="intersects") contacts = contacts[["X", "Y", "Z", "geometry", "basal_unit"]].copy() - bounding_box = map_data.get_bounding_box() - # Interpolate the dip of the contacts interpolator = DipDipDirectionInterpolator(data_type="dip") # Interpolate the dip of the contacts - dip = interpolator(bounding_box, structure_data, interpolator=scipy.interpolate.Rbf) + dip = interpolator(self.bounding_box, structure_data, interpolator=scipy.interpolate.Rbf) # create a GeoDataFrame of the interpolated orientations interpolated_orientations = geopandas.GeoDataFrame() # add the dip and dip direction to the GeoDataFrame @@ -301,14 +316,13 @@ def compute( # set the crs of the interpolated orientations to the crs of the units interpolated_orientations = interpolated_orientations.set_crs(crs=basal_contacts.crs) # get the elevation Z of the interpolated points - dtm_data = map_data.get_map_data(Datatype.DTM) - interpolated = set_z_values_from_raster_df(dtm_data, interpolated_orientations) + interpolated = set_z_values_from_raster_df(self.dtm_data, interpolated_orientations) # update the geometry of the interpolated points to include the Z value interpolated["geometry"] = interpolated.apply( lambda row: shapely.geometry.Point(row.geometry.x, row.geometry.y, row["Z"]), axis=1 ) # for each interpolated point, assign name of unit using spatial join - units = map_data.get_map_data(Datatype.GEOLOGY) + units = geology_data.copy() interpolated_orientations = interpolated_orientations.sjoin( units, how="inner", predicate="within" ) @@ -358,13 +372,13 @@ def compute( p1[0] = numpy.asarray(short_line[0].coords[0][0]) p1[1] = numpy.asarray(short_line[0].coords[0][1]) # get the elevation Z of the end point p1 - p1[2] = map_data.get_value_from_raster(Datatype.DTM, p1[0], p1[1]) + p1[2] = value_from_raster(Datatype.DTM, p1[0], p1[1]) # create array to store xyz coordinates of the end point p2 p2 = numpy.zeros(3) p2[0] = numpy.asarray(short_line[0].coords[-1][0]) p2[1] = numpy.asarray(short_line[0].coords[-1][1]) # get the elevation Z of the end point p2 - p2[2] = map_data.get_value_from_raster(Datatype.DTM, p2[0], p2[1]) + p2[2] = value_from_raster(Datatype.DTM, p2[0], p2[1]) # calculate the length of the shortest line line_length = scipy.spatial.distance.euclidean(p1, p2) # find the indices of the points that are within 5% of the length of the shortest line @@ -444,8 +458,13 @@ class StructuralPoint(ThicknessCalculator): ''' - def __init__(self, max_line_length: float = None): - super().__init__(max_line_length) + def __init__( + self, + dtm_data: Optional[gdal.Dataset] = None, + bounding_box: Optional[dict] = None, + max_line_length: Optional[float] = None + ): + super().__init__(dtm_data, bounding_box, max_line_length) self.thickness_calculator_label = "StructuralPoint" self.strike_allowance = 30 self.lines = None @@ -458,7 +477,8 @@ def compute( stratigraphic_order: list, basal_contacts: geopandas.GeoDataFrame, structure_data: pandas.DataFrame, - map_data: MapData, + geology_data: geopandas.GeoDataFrame, + sampled_contacts: geopandas.GeoDataFrame, ) -> pandas.DataFrame: """ Method overview: @@ -499,7 +519,7 @@ def compute( basal_contacts = basal_contacts.copy() # grab geology polygons and calculate bounding boxes for each lithology - geology = map_data.get_map_data(datatype=Datatype.GEOLOGY) + geology = geology_data.copy() geology[['minx', 'miny', 'maxx', 'maxy']] = geology.bounds # create a GeoDataFrame of the sampled structures @@ -520,7 +540,7 @@ def compute( # rebuild basal contacts lines based on sampled dataset sampled_basal_contacts = rebuild_sampled_basal_contacts( - basal_contacts, map_data.sampled_contacts + basal_contacts, sampled_contacts ) # calculate map dimensions From 4bfdbb07b90edf99314b84863f118b9ae3ba866c Mon Sep 17 00:00:00 2001 From: rabii-chaarani Date: Mon, 28 Jul 2025 13:17:58 +0930 Subject: [PATCH 2/5] refactor: update initialisation of thickness calculator --- map2loop/project.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/map2loop/project.py b/map2loop/project.py index 0dc62b0b..ab053085 100644 --- a/map2loop/project.py +++ b/map2loop/project.py @@ -41,7 +41,7 @@ class Project(object): verbose_level: m2l_enums.VerboseLevel A selection that defines how much console logging is output samplers: Sampler - A list of samplers used to extract point samples from polyonal or line segments. Indexed by m2l_enum.Dataype + A list of samplers used to extract point samples from polygons or line segments. Indexed by m2l_enum sorter: Sorter The sorting algorithm to use for calculating the stratigraphic column loop_filename: str @@ -141,7 +141,6 @@ def __init__( self.bounding_box = bounding_box self.contact_extractor = None self.sorter = SorterUseHint() - self.thickness_calculator = [InterpolatedStructure()] self.throw_calculator = ThrowCalculatorAlpha() self.fault_orientation = FaultOrientationNearest() self.map_data = MapData(verbose_level=verbose_level) @@ -245,6 +244,10 @@ def __init__( self.map_data.get_map_data(Datatype.GEOLOGY), self.map_data.get_map_data(Datatype.FAULT) ) + self.thickness_calculator = [InterpolatedStructure( + dtm_data=self.map_data.get_map_data(Datatype.DTM), + bounding_box=self.bounding_box, + )] @beartype.beartype @@ -697,7 +700,8 @@ def calculate_unit_thicknesses(self, basal_contacts): self.stratigraphic_column.column, basal_contacts, self.structure_samples, - self.map_data, + self.map_data.get_map_data(Datatype.GEOLOGY), + self.map_data.sampled_contacts, )[['ThicknessMean', 'ThicknessMedian', 'ThicknessStdDev']].to_numpy() label = calculator.thickness_calculator_label From 7ec4d314f5ec9685563a5e9379be2adcc2d4cfe9 Mon Sep 17 00:00:00 2001 From: rabii-chaarani Date: Mon, 28 Jul 2025 15:10:27 +0930 Subject: [PATCH 3/5] refactor: update tests --- .../InterpolatedStructure/test_interpolated_structure.py | 6 ++++-- .../StructurePoint/test_ThicknessStructuralPoint.py | 3 ++- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/tests/thickness/InterpolatedStructure/test_interpolated_structure.py b/tests/thickness/InterpolatedStructure/test_interpolated_structure.py index 9bc28431..12d0c87c 100644 --- a/tests/thickness/InterpolatedStructure/test_interpolated_structure.py +++ b/tests/thickness/InterpolatedStructure/test_interpolated_structure.py @@ -1663,7 +1663,6 @@ def check_thickness_values(result, column, description): def test_calculate_thickness_InterpolatedStructure(): # Run the calculation - thickness_calculator = InterpolatedStructure() md = MapData() md.sampled_contacts = s_c @@ -1683,13 +1682,16 @@ def test_calculate_thickness_InterpolatedStructure(): "base": -3200, "top": 3000, } + thickness_calculator = InterpolatedStructure(dtm_data=md.get_map_data(Datatype.DTM), + bounding_box=md.bounding_box) result = thickness_calculator.compute( units=st_units, stratigraphic_order=st_column, basal_contacts=bc_gdf, structure_data=structures, - map_data=md, + geology_data=md.get_map_data(Datatype.GEOLOGY), + sampled_contacts=md.sampled_contacts, ) # is thickness calc alpha the label? diff --git a/tests/thickness/StructurePoint/test_ThicknessStructuralPoint.py b/tests/thickness/StructurePoint/test_ThicknessStructuralPoint.py index a71ebb7b..b7e2e7d8 100644 --- a/tests/thickness/StructurePoint/test_ThicknessStructuralPoint.py +++ b/tests/thickness/StructurePoint/test_ThicknessStructuralPoint.py @@ -1675,7 +1675,8 @@ def test_calculate_thickness_structural_point(): stratigraphic_order=st_column, basal_contacts=bc_gdf, structure_data=structures, - map_data=md, + geology_data=md.get_map_data(Datatype.GEOLOGY), + sampled_contacts=md.sampled_contacts, ) # is thickness calc alpha the label? From 2a678e53b015897796620551c812a010cdd3a926 Mon Sep 17 00:00:00 2001 From: Noelle Cheng Date: Fri, 1 Aug 2025 13:07:03 +0800 Subject: [PATCH 4/5] fix parameters pass to value_from_raster in thickness calculator --- map2loop/thickness_calculator.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/map2loop/thickness_calculator.py b/map2loop/thickness_calculator.py index f074f9b3..b952a6c4 100644 --- a/map2loop/thickness_calculator.py +++ b/map2loop/thickness_calculator.py @@ -366,19 +366,22 @@ def compute( # check if the short line is if self.max_line_length is not None and short_line.length > self.max_line_length: continue + + inv_geotransform = gdal.InvGeoTransform(self.dtm_data.GetGeoTransform()) + data_array = numpy.array(self.dtm_data.GetRasterBand(1).ReadAsArray().T) # extract the end points of the shortest line p1 = numpy.zeros(3) p1[0] = numpy.asarray(short_line[0].coords[0][0]) p1[1] = numpy.asarray(short_line[0].coords[0][1]) # get the elevation Z of the end point p1 - p1[2] = value_from_raster(Datatype.DTM, p1[0], p1[1]) + p1[2] = value_from_raster(inv_geotransform, data_array, p1[0], p1[1]) # create array to store xyz coordinates of the end point p2 p2 = numpy.zeros(3) p2[0] = numpy.asarray(short_line[0].coords[-1][0]) p2[1] = numpy.asarray(short_line[0].coords[-1][1]) # get the elevation Z of the end point p2 - p2[2] = value_from_raster(Datatype.DTM, p2[0], p2[1]) + p2[2] = value_from_raster(inv_geotransform, data_array, p2[0], p2[1]) # calculate the length of the shortest line line_length = scipy.spatial.distance.euclidean(p1, p2) # find the indices of the points that are within 5% of the length of the shortest line From 5193b5a3e76059ecdf2c4c6f17cd9dd4c237dafc Mon Sep 17 00:00:00 2001 From: Noelle Cheng Date: Fri, 1 Aug 2025 13:08:13 +0800 Subject: [PATCH 5/5] fix thickness calculator tests --- .../test_interpolated_structure.py | 4 +++- .../test_ThicknessStructuralPoint.py | 5 +++-- .../test_ThicknessCalculatorAlpha.py | 17 ++++++++++------- 3 files changed, 16 insertions(+), 10 deletions(-) diff --git a/tests/thickness/InterpolatedStructure/test_interpolated_structure.py b/tests/thickness/InterpolatedStructure/test_interpolated_structure.py index 12d0c87c..1a0f65c4 100644 --- a/tests/thickness/InterpolatedStructure/test_interpolated_structure.py +++ b/tests/thickness/InterpolatedStructure/test_interpolated_structure.py @@ -1,6 +1,7 @@ import pandas import geopandas import numpy +import shapely.geometry from map2loop.mapdata import MapData from map2loop.thickness_calculator import InterpolatedStructure @@ -1640,7 +1641,8 @@ '3', ] -s_c = pandas.DataFrame({'X': X, 'Y': Y, 'Z': Z, 'featureId': featureid}) +geometry= [shapely.geometry.Point(x,y) for x,y in zip(X, Y)] +s_c = geopandas.GeoDataFrame({'X': X, 'Y': Y, 'Z': Z, 'featureId': featureid}, geometry = geometry, crs='EPSG:28350') ################################## diff --git a/tests/thickness/StructurePoint/test_ThicknessStructuralPoint.py b/tests/thickness/StructurePoint/test_ThicknessStructuralPoint.py index b7e2e7d8..08d0aabf 100644 --- a/tests/thickness/StructurePoint/test_ThicknessStructuralPoint.py +++ b/tests/thickness/StructurePoint/test_ThicknessStructuralPoint.py @@ -1,6 +1,7 @@ import pandas import geopandas import numpy +import shapely.geometry from map2loop.mapdata import MapData from map2loop.thickness_calculator import StructuralPoint @@ -1637,7 +1638,8 @@ '3', ] -s_c = pandas.DataFrame({'X': X, 'Y': Y, 'Z': Z, 'featureId': featureid}) +geometry= [shapely.geometry.Point(x,y) for x,y in zip(X, Y)] +s_c = geopandas.GeoDataFrame({'X': X, 'Y': Y, 'Z': Z, 'featureId': featureid}, geometry = geometry, crs='EPSG:28350') ############################ @@ -1664,7 +1666,6 @@ def test_calculate_thickness_structural_point(): md = MapData() md.sampled_contacts = s_c - md.sampled_contacts = s_c md.raw_data[Datatype.GEOLOGY] = geology md.load_map_data(Datatype.GEOLOGY) md.check_map(Datatype.GEOLOGY) diff --git a/tests/thickness/ThicknessCalculatorAlpha/test_ThicknessCalculatorAlpha.py b/tests/thickness/ThicknessCalculatorAlpha/test_ThicknessCalculatorAlpha.py index a169b3ce..827a1859 100644 --- a/tests/thickness/ThicknessCalculatorAlpha/test_ThicknessCalculatorAlpha.py +++ b/tests/thickness/ThicknessCalculatorAlpha/test_ThicknessCalculatorAlpha.py @@ -1,6 +1,7 @@ import pandas import geopandas import numpy +import shapely.geometry from map2loop.mapdata import MapData from map2loop.m2l_enums import Datatype @@ -1639,8 +1640,8 @@ '3', ] -s_c = pandas.DataFrame({'X': X, 'Y': Y, 'Z': Z, 'featureId': featureid}) - +geometry= [shapely.geometry.Point(x,y) for x,y in zip(X, Y)] +s_c = geopandas.GeoDataFrame({'X': X, 'Y': Y, 'Z': Z, 'featureId': featureid}, geometry = geometry, crs='EPSG:28350') ##################################### ### TEST ThicknessCalculatorAlpha ### @@ -1658,15 +1659,16 @@ def check_thickness_values(result, column, description): geology = load_hamersley_geology() - +geology.rename(columns={'unitname': 'UNITNAME', 'code': 'CODE'}, inplace=True) def test_calculate_thickness_thickness_calculator_alpha(): # Run the calculation md = MapData() md.sampled_contacts = s_c - md.data[Datatype.GEOLOGY] = geology - - print('GERE', md.get_map_data(Datatype.GEOLOGY)) + md.raw_data[Datatype.GEOLOGY] = geology + md.load_map_data(Datatype.GEOLOGY) + md.check_map(Datatype.GEOLOGY) + md.parse_geology_map() thickness_calculator = ThicknessCalculatorAlpha() result = thickness_calculator.compute( @@ -1674,7 +1676,8 @@ def test_calculate_thickness_thickness_calculator_alpha(): stratigraphic_order=st_column, basal_contacts=bc_gdf, structure_data=structures, - map_data=md, + geology_data=md.get_map_data(Datatype.GEOLOGY), + sampled_contacts=md.sampled_contacts, ) # is thickness calc alpha the label?