Thickness Calculators - Removed Mapdata Dependency

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

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,29 +282,26 @@ 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)
 
         # 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
         )
         # 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"
         )
@@ -352,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] = map_data.get_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] = map_data.get_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
@@ -444,8 +461,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 +480,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 +522,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 +543,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

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')
 
 
 ##################################
@@ -1663,7 +1665,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 +1684,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?

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)
@@ -1675,7 +1676,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?

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,23 +1659,25 @@ 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(
         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?