From d936aa592441bb698849e7975f300d01e6973c68 Mon Sep 17 00:00:00 2001
From: rabii-chaarani <rabii.chaarani@monash.edu>
Date: Mon, 10 Nov 2025 12:49:44 +0930
Subject: [PATCH 01/31] feat: add AlongSection thickness calculator

---
 map2loop/thickness_calculator.py | 314 ++++++++++++++++++++++++++++++-
 1 file changed, 312 insertions(+), 2 deletions(-)

diff --git a/map2loop/thickness_calculator.py b/map2loop/thickness_calculator.py
index 13a1b806..0e48547d 100644
--- a/map2loop/thickness_calculator.py
+++ b/map2loop/thickness_calculator.py
@@ -6,7 +6,11 @@
     multiline_to_line,
     find_segment_strike_from_pt,
     set_z_values_from_raster_df,
-    value_from_raster
+    value_from_raster,
+    segment_measure_range,
+    clean_line_geometry,
+    nearest_orientation_to_line,
+    iter_line_segments
 )
 from .interpolators import DipDipDirectionInterpolator
 
@@ -15,11 +19,13 @@
 
 # external imports
 from abc import ABC, abstractmethod
-from typing import Optional
+from typing import Optional, List
+from collections import defaultdict
 import scipy.interpolate
 import beartype
 import numpy
 import scipy
+from scipy.spatial import cKDTree
 import pandas
 import geopandas
 import shapely
@@ -756,3 +762,307 @@ def compute(
         self._check_thickness_percentage_calculations(output_units)
         
         return output_units
+
+
+class AlongSection(ThicknessCalculator):
+    """Thickness calculator that estimates true thicknesses along supplied section lines."""
+
+    def __init__(
+        self,
+        sections: geopandas.GeoDataFrame,
+        dtm_data: Optional[gdal.Dataset] = None,
+        bounding_box: Optional[dict] = None,
+        max_line_length: Optional[float] = None,
+        is_strike: Optional[bool] = False,
+    ):
+        super().__init__(dtm_data, bounding_box, max_line_length, is_strike)  # initialise base calculator bits
+        self.thickness_calculator_label = "AlongSection"  # label used externally to identify this strategy
+        self.sections = (
+            sections.copy()  # keep a copy so editing outside does not mutate our internal state
+            if sections is not None
+            else geopandas.GeoDataFrame({"geometry": []}, geometry="geometry")  # ensure predictable structure when no sections supplied
+        )
+        self.section_thickness_records: geopandas.GeoDataFrame | List[dict] = []  # populated as GeoDataFrame during compute
+        self.section_intersection_points: dict = {}
+
+    def type(self):
+        """Return the calculator label."""
+        return self.thickness_calculator_label
+
+    @beartype.beartype
+    def compute(
+        self,
+        units: pandas.DataFrame,
+        stratigraphic_order: list,
+        basal_contacts: geopandas.GeoDataFrame,
+        structure_data: pandas.DataFrame,
+        geology_data: geopandas.GeoDataFrame,
+        sampled_contacts: pandas.DataFrame,
+    ) -> pandas.DataFrame:
+        """Estimate unit thicknesses along cross sections.
+
+        Args:
+            units (pandas.DataFrame): Table of stratigraphic units that will be annotated with
+                thickness statistics. Must expose a ``name`` column used to map aggregated
+                thickness values back to each unit.
+            stratigraphic_order (list): Accepted for API compatibility but not used by this
+                calculator. Retained so the method signature matches the other calculators.
+            basal_contacts (geopandas.GeoDataFrame): Unused placeholder maintained for parity
+                with the other calculators.
+            structure_data (pandas.DataFrame): Orientation measurements containing X/Y columns
+                and one dip column (``DIP``, ``dip`` or ``Dip``). A KD-tree is built from these
+                points so each split section segment can grab the nearest dip when converting its
+                length to true thickness.
+            geology_data (geopandas.GeoDataFrame): Mandatory polygon dataset describing map
+                units. Each section is intersected with these polygons to generate unit-aware
+                line segments whose lengths underpin the thickness calculation.
+            sampled_contacts (pandas.DataFrame): Part of the public interface but not used in
+                this implementation.
+
+        Returns:
+            pandas.DataFrame: Copy of ``units`` with three new columns – ``ThicknessMean``,
+            ``ThicknessMedian`` and ``ThicknessStdDev`` – populated from the accumulated segment
+            thickness samples. Units that never receive a measurement retain the sentinel value
+            ``-1`` in each column.
+
+        Workflow:
+            1. Validate the presence of section geometries, geology polygons, and a recognizable
+               unit-name column. Reproject sections to match the geology CRS when required.
+            2. Pre-build helpers: a spatial index over geology polygons for fast section/polygon
+               queries, and (if possible) a KD-tree of orientation points so the closest dip to a
+               split segment can be queried in O(log n).
+            3. For every section, reduce complex/multi-part geometries to a single `LineString`,
+               intersect it with candidate geology polygons, and collect the resulting split line
+               segments along with their length and measure positions along the parent section.
+            4. Walk the ordered segments and keep those bounded by two distinct neighbouring units
+               (i.e., the segment sits between different units on either side). Fetch the nearest
+               dip (fallback to 90° once if no structures exist), convert the segment length to
+               true thickness using ``length * sin(dip)``, and store the sample plus provenance in
+               ``self.section_thickness_records``.
+            5. Aggregate all collected samples per unit (mean/median/std) and return the enriched
+               ``units`` table. The helper ``self.section_intersection_points`` captures the raw
+               split segments grouped by section to aid downstream inspection.
+
+        Notes:
+            - ``stratigraphic_order``, ``basal_contacts`` and ``sampled_contacts`` are retained for
+              interface compatibility but do not influence this algorithm.
+            - When no nearby orientation measurements exist, the method emits a single warning and
+              assumes a vertical dip (90°) for affected segments to avoid silently dropping units.
+            - ``self.section_thickness_records`` is materialised as a ``GeoDataFrame`` so the split
+              segments (geometry column) can be visualised directly in notebooks or GIS clients.
+        """
+
+        if self.sections is None or self.sections.empty:
+            logger.warning("AlongSection: No sections provided; skipping thickness calculation.")
+            return units
+
+        if geology_data is None or geology_data.empty:
+            logger.warning(
+                "AlongSection: Geology polygons are required to split sections; skipping thickness calculation."
+            )
+            return units
+
+        unit_column_candidates = [
+            "UNITNAME",
+            "unitname",
+            "UNIT_NAME",
+            "UnitName",
+            "unit",
+            "Unit",
+            "UNIT",
+            "name",
+            "Name",
+        ]
+        unit_column = next((col for col in unit_column_candidates if col in geology_data.columns), None)
+        if unit_column is None:
+            logger.warning(
+                "AlongSection: Unable to identify a unit-name column in geology data; expected one of %s.",
+                unit_column_candidates,
+            )
+            return units
+
+        sections = self.sections.copy()
+        geology = geology_data.copy()
+
+        if geology.crs is not None:
+            if sections.crs is None:
+                sections = sections.set_crs(geology.crs)
+            elif sections.crs != geology.crs:
+                sections = sections.to_crs(geology.crs)
+        elif sections.crs is not None:
+            geology = geology.set_crs(sections.crs)
+
+        try:
+            geology_sindex = geology.sindex
+        except Exception:
+            geology_sindex = None
+
+        dip_column = next((col for col in ("DIP", "dip", "Dip") if col in structure_data.columns), None)
+        orientation_tree = None
+        orientation_dips = None
+        orientation_coords = None
+        if dip_column is not None and {"X", "Y"}.issubset(structure_data.columns):
+            orient_df = structure_data.dropna(subset=["X", "Y", dip_column]).copy()
+            if not orient_df.empty:
+                orient_df["_dip_value"] = pandas.to_numeric(orient_df[dip_column], errors="coerce")
+                orient_df = orient_df.dropna(subset=["_dip_value"])
+                if not orient_df.empty:
+                    try:
+                        orientation_coords = orient_df[["X", "Y"]].astype(float).to_numpy()
+                    except ValueError:
+                        orientation_coords = numpy.empty((0, 2))
+                    if orientation_coords.size:
+                        orientation_dips = orient_df["_dip_value"].astype(float).to_numpy()
+                        try:
+                            orientation_tree = cKDTree(orientation_coords)
+                        except Exception:
+                            orientation_tree = None
+
+        default_dip_warning_emitted = False
+
+        units_lookup = dict(zip(units["name"], units.index))
+        thickness_by_unit = {name: [] for name in units_lookup.keys()}
+        thickness_records: List[dict] = []
+        split_segments_by_section: dict = defaultdict(list)
+
+        for section_idx, section_row in sections.iterrows():
+            line = clean_line_geometry(section_row.geometry)
+            if line is None or line.length == 0:
+                continue
+            section_id = section_row.get("ID", section_idx)
+
+            candidate_idx = (
+                list(geology_sindex.intersection(line.bounds))
+                if geology_sindex is not None
+                else list(range(len(geology)))
+            )
+            if not candidate_idx:
+                continue
+
+            split_segments = []
+            for _, poly in geology.iloc[candidate_idx].iterrows():
+                polygon_geom = poly.geometry
+                if polygon_geom is None or polygon_geom.is_empty:
+                    continue
+                intersection = line.intersection(polygon_geom)
+                if intersection.is_empty:
+                    continue
+                for segment in iter_line_segments(intersection):
+                    seg_length = segment.length
+                    if seg_length <= 0:
+                        continue
+                    start_measure, end_measure = segment_measure_range(line, segment)
+                    segment_record = {
+                        "section_id": section_id,
+                        "geometry": segment,
+                        "unit": poly[unit_column],
+                        "length": seg_length,
+                        "start_measure": start_measure,
+                        "end_measure": end_measure,
+                    }
+                    split_segments.append(segment_record)
+                    split_segments_by_section[section_id].append(segment_record)
+
+            if len(split_segments) < 1:
+                continue
+
+            split_segments.sort(key=lambda item: (item["start_measure"], item["end_measure"]))
+
+            for idx_segment, segment in enumerate(split_segments):
+                prev_unit = split_segments[idx_segment - 1]["unit"] if idx_segment > 0 else None
+                next_unit = (
+                    split_segments[idx_segment + 1]["unit"]
+                    if idx_segment + 1 < len(split_segments)
+                    else None
+                )
+
+                if prev_unit is None or next_unit is None:
+                    continue
+                if prev_unit == segment["unit"] or next_unit == segment["unit"]:
+                    continue
+                if prev_unit == next_unit:
+                    continue
+
+                dip_value, orientation_idx, orientation_distance = nearest_orientation_to_line(
+                    orientation_tree, 
+                    orientation_dips,
+                    orientation_coords,
+                    segment["geometry"]
+                )
+                if numpy.isnan(dip_value):
+                    dip_value = 90.0
+                    orientation_idx = None
+                    orientation_distance = None
+                    if not default_dip_warning_emitted:
+                        logger.warning(
+                            "AlongSection: Missing structure measurements near some sections; assuming vertical dip (90°) for those segments."
+                        )
+                        default_dip_warning_emitted = True
+
+                thickness = segment["length"] * abs(math.sin(math.radians(dip_value)))
+                if thickness <= 0:
+                    continue
+
+                unit_name = segment["unit"]
+                if unit_name not in thickness_by_unit:
+                    thickness_by_unit[unit_name] = []
+                thickness_by_unit[unit_name].append(thickness)
+
+                thickness_records.append(
+                    {
+                        "section_id": section_id,
+                        "unit_id": unit_name,
+                        "thickness": thickness,
+                        "segment_length": segment["length"],
+                        "dip_used_deg": dip_value,
+                        "prev_unit": prev_unit,
+                        "next_unit": next_unit,
+                        "orientation_index": orientation_idx,
+                        "orientation_distance": orientation_distance,
+                        "geometry": segment["geometry"],
+                    }
+                )
+
+        output_units = units.copy()
+        output_units["ThicknessMean"] = -1.0
+        output_units["ThicknessMedian"] = -1.0
+        output_units["ThicknessStdDev"] = -1.0
+
+        for unit_name, values in thickness_by_unit.items():
+            if not values:
+                continue
+            idx = units_lookup.get(unit_name)
+            if idx is None:
+                continue
+            arr = numpy.asarray(values, dtype=numpy.float64)
+            output_units.at[idx, "ThicknessMean"] = float(numpy.nanmean(arr))
+            output_units.at[idx, "ThicknessMedian"] = float(numpy.nanmedian(arr))
+            output_units.at[idx, "ThicknessStdDev"] = float(numpy.nanstd(arr))
+
+        if thickness_records:
+            self.section_thickness_records = geopandas.GeoDataFrame(
+                thickness_records,
+                geometry="geometry",
+                crs=sections.crs,
+            )
+        else:
+            self.section_thickness_records = geopandas.GeoDataFrame(
+                columns=[
+                    "section_id",
+                    "unit_id",
+                    "thickness",
+                    "segment_length",
+                    "dip_used_deg",
+                    "prev_unit",
+                    "next_unit",
+                    "orientation_index",
+                    "orientation_distance",
+                    "geometry",
+                ],
+                geometry="geometry",
+                crs=sections.crs,
+            )
+        self.section_intersection_points = dict(split_segments_by_section)
+        self._check_thickness_percentage_calculations(output_units)
+
+        return output_units

From 3822263b64bf646e07c0808042dbb14ff95e4664 Mon Sep 17 00:00:00 2001
From: rabii-chaarani <rabii.chaarani@monash.edu>
Date: Mon, 10 Nov 2025 12:50:17 +0930
Subject: [PATCH 02/31] feat: add utils for AlongSection class

---
 map2loop/utils.py | 209 +++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 208 insertions(+), 1 deletion(-)

diff --git a/map2loop/utils.py b/map2loop/utils.py
index 55e2e7b2..ce84ee60 100644
--- a/map2loop/utils.py
+++ b/map2loop/utils.py
@@ -587,4 +587,211 @@ def set_z_values_from_raster_df(dtm_data: gdal.Dataset, df: pandas.DataFrame):
         axis=1,
     )
 
-    return df
\ No newline at end of file
+    return df
+
+def clean_line_geometry(geometry: shapely.geometry.base.BaseGeometry):
+    """Clean and normalize a Shapely geometry to a single LineString or None.
+
+    Args:
+        geometry (shapely.geometry.base.BaseGeometry | None): Input geometry to be
+            cleaned. Accepted inputs include ``LineString``, ``MultiLineString``,
+            ``GeometryCollection``, and other Shapely geometries. If ``None`` or
+            an empty geometry is provided, the function returns ``None``.
+
+    Returns:
+        shapely.geometry.LineString or None: A single ``LineString`` if the
+        geometry is already a ``LineString`` or can be merged/converted into one.
+        If the input is ``None``, empty, or cannot be converted into a valid
+        ``LineString``, the function returns ``None``.
+
+    Raises:
+        None: Exceptions raised by ``shapely.ops.linemerge`` are caught and result
+        in a ``None`` return rather than being propagated.
+
+    Notes:
+        - The function uses ``shapely.ops.linemerge`` to attempt to merge multipart
+          geometries into a single ``LineString``.
+        - If ``linemerge`` returns a ``MultiLineString``, the longest constituent
+          ``LineString`` (by ``length``) is returned.
+        - Geometries that are already ``LineString`` are returned unchanged.
+
+    Examples:
+        >>> from shapely.geometry import LineString
+        >>> clean_line_geometry(None) is None
+        True
+        >>> ls = LineString([(0, 0), (1, 1)])
+        >>> clean_line_geometry(ls) is ls
+        True
+    """
+    if geometry is None or geometry.is_empty:
+        return None
+    if geometry.geom_type == "LineString":
+        return geometry
+    try:
+        merged = shapely.ops.linemerge(geometry)
+    except Exception:
+        return None
+    if merged.geom_type == "LineString":
+        return merged
+    if merged.geom_type == "MultiLineString":
+        try:
+            return max(merged.geoms, key=lambda geom: geom.length)
+        except ValueError:
+            return None
+    return None
+
+def iter_line_segments(geometry):
+    """Produce all LineString segments contained in a Shapely geometry.
+
+    Args:
+        geometry (shapely.geometry.base.BaseGeometry | None): Input geometry. Accepted
+            types include ``LineString``, ``MultiLineString``, ``GeometryCollection``
+            and other Shapely geometries that may contain line parts. If ``None`` or
+            an empty geometry is provided, an empty list is returned.
+
+    Returns:
+        list[shapely.geometry.LineString]: A list of ``LineString`` objects extracted
+        from the input geometry. Behavior by input type:
+            - ``LineString``: returned as a single-element list.
+            - ``MultiLineString``: returns the non-zero-length constituent parts.
+            - ``GeometryCollection``: recursively extracts contained line segments.
+            - Other geometry types: returns an empty list if no line segments found.
+
+    Notes:
+        Zero-length segments are filtered out. The function is defensive and
+        will return an empty list for ``None`` or empty geometries rather than
+        raising an exception.
+
+    Examples:
+        >>> from shapely.geometry import LineString
+        >>> iter_line_segments(LineString([(0, 0), (1, 1)]))
+        [LineString([(0, 0), (1, 1)])]
+    """
+    if geometry is None or geometry.is_empty:
+        return []
+    if isinstance(geometry, shapely.geometry.LineString):
+        return [geometry]
+    if isinstance(geometry, shapely.geometry.MultiLineString):
+        return [geom for geom in geometry.geoms if geom.length > 0]
+    if isinstance(geometry, shapely.geometry.GeometryCollection):
+        segments = []
+        for geom in geometry.geoms:
+            segments.extend(iter_line_segments(geom))
+        return segments
+    return []
+
+def nearest_orientation_to_line(orientation_tree, orientation_dips, orientation_coords, line_geom: shapely.geometry.LineString):
+    """Find the nearest orientation measurement to the midpoint of a line.
+
+    This function queries the provided spatial index and orientation arrays and
+    returns the dip value, the index of the matched orientation, and the
+    distance from the line to that orientation point.
+
+    Args:
+        orientation_tree: Spatial index object (e.g., KDTree) with a ``query``
+            method accepting an (x, y) tuple and optional ``k``. Typically
+            provided by the caller so the utility remains stateless.
+        orientation_dips (Sequence[float]): Sequence of dip values aligned with
+            ``orientation_coords``.
+        orientation_coords (Sequence[tuple]): Sequence of (x, y) coordinate
+            tuples for each orientation measurement.
+        line_geom (shapely.geometry.LineString): Line geometry for which the
+            nearest orientation should be found. The midpoint (interpolated at
+            50% of the line length) is used as the query point; if the
+            midpoint is empty, the line centroid is used as a fallback.
+
+    Returns:
+        tuple: A tuple ``(dip, index, distance)`` where:
+            - ``dip`` (float or numpy.nan): the dip value from
+              ``orientation_dips`` at the nearest orientation. Returns
+              ``numpy.nan`` if no valid orientation can be found or on error.
+            - ``index`` (int or None): the integer index into the orientation
+              arrays for the best match, or ``None`` if not found.
+            - ``distance`` (float or None): the shortest geometric distance
+              between ``line_geom`` and the matched orientation point, or
+              ``None`` if not available.
+
+    Notes:
+        - The function queries up to 5 nearest neighbors (or fewer if fewer
+          orientations exist) and then computes the actual geometry distance
+          between the ``line_geom`` and each candidate point to select the
+          closest match. Any exceptions from the spatial query are caught and
+          result in ``(numpy.nan, None, None)`` being returned instead of
+          propagating the exception.
+        - The function is defensive: invalid or empty geometries return
+          ``(numpy.nan, None, None)`` rather than raising.
+
+    Examples:
+        >>> dip, idx, dist = nearest_orientation_to_line(tree, dips, coords, my_line)
+        >>> if idx is not None:
+        ...     print(f"Nearest dip={dip} at index={idx} (distance={dist})")
+    """
+    if orientation_tree is None or orientation_dips is None or orientation_coords is None:
+        return numpy.nan, None, None
+    midpoint = line_geom.interpolate(0.5, normalized=True)
+    if midpoint.is_empty:
+        midpoint = line_geom.centroid
+    if midpoint.is_empty:
+        return numpy.nan, None, None
+    query_xy = (float(midpoint.x), float(midpoint.y))
+    k = min(len(orientation_dips), 5)
+    try:
+        distances, indices = orientation_tree.query(query_xy, k=k)
+    except Exception:
+        return numpy.nan, None, None
+    distances = numpy.atleast_1d(distances)
+    indices = numpy.atleast_1d(indices)
+    best_idx = None
+    best_dist = numpy.inf
+    for approx_dist, idx in zip(distances, indices):
+        if idx is None:
+            continue
+        candidate_point = shapely.geometry.Point(orientation_coords[int(idx)])
+        actual_dist = line_geom.distance(candidate_point)
+        if actual_dist < best_dist:
+            best_dist = actual_dist
+            best_idx = int(idx)
+    if best_idx is None:
+        return numpy.nan, None, None
+    return float(orientation_dips[best_idx]), best_idx, best_dist
+
+def segment_measure_range(parent_line: shapely.geometry.LineString, segment: shapely.geometry.LineString):
+    """Compute projected measures of a segment's end points along a parent line.
+
+    This function projects the start and end points of ``segment`` onto
+    ``parent_line`` using Shapely's ``project`` method and returns the two
+    measures in ascending order (start <= end). Measures are distances along
+    the parent line's linear reference (units of the geometry's CRS).
+
+    Args:
+        parent_line (shapely.geometry.LineString): The reference line onto which
+            the segment end points will be projected.
+        segment (shapely.geometry.LineString): The segment whose first and last
+            vertices will be projected onto ``parent_line``. The function uses
+            ``segment.coords[0]`` and ``segment.coords[-1]`` as the start and
+            end points respectively.
+
+    Returns:
+        tuple(float, float): A tuple ``(start_measure, end_measure)`` containing
+        the projected distances along ``parent_line`` for the segment's start
+        and end points. The values are ordered so that ``start_measure <= end_measure``.
+
+    Notes:
+        - If ``segment`` is degenerate (e.g., a single point), both measures may
+          be equal.
+        - The returned measures are in the same linear units as the input
+          geometries (e.g., metres for projected CRS).
+        - No validation is performed on the inputs; callers should ensure both
+          geometries are valid and share a common CRS.
+
+    Examples:
+        >>> start_m, end_m = segment_measure_range(parent_line, segment)
+        >>> assert start_m <= end_m
+    """
+    start_point = shapely.geometry.Point(segment.coords[0])
+    end_point = shapely.geometry.Point(segment.coords[-1])
+    start_measure = parent_line.project(start_point)
+    end_measure = parent_line.project(end_point)
+    if end_measure < start_measure:
+        start_measure, end_measure = end_measure, start_measure
+    return start_measure, end_measure
\ No newline at end of file

From 7a2223405833a4c7f398ba92ec52134a21084970 Mon Sep 17 00:00:00 2001
From: rabii-chaarani <rabii.chaarani@monash.edu>
Date: Mon, 10 Nov 2025 12:54:22 +0930
Subject: [PATCH 03/31] fix: renamed to _approx_dist

---
 map2loop/utils.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/map2loop/utils.py b/map2loop/utils.py
index ce84ee60..927d4c5e 100644
--- a/map2loop/utils.py
+++ b/map2loop/utils.py
@@ -743,7 +743,7 @@ def nearest_orientation_to_line(orientation_tree, orientation_dips, orientation_
     indices = numpy.atleast_1d(indices)
     best_idx = None
     best_dist = numpy.inf
-    for approx_dist, idx in zip(distances, indices):
+    for _approx_dist, idx in zip(distances, indices):
         if idx is None:
             continue
         candidate_point = shapely.geometry.Point(orientation_coords[int(idx)])

From aa813aa13fc2bbd3b9b4603d52e921f0b4372d4b Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Thu, 13 Nov 2025 12:35:19 +1100
Subject: [PATCH 04/31] Adding exception handling suggestions

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
---
 map2loop/thickness_calculator.py | 17 ++++++++++++-----
 map2loop/utils.py                |  6 ++++--
 2 files changed, 16 insertions(+), 7 deletions(-)

diff --git a/map2loop/thickness_calculator.py b/map2loop/thickness_calculator.py
index 0e48547d..614344cb 100644
--- a/map2loop/thickness_calculator.py
+++ b/map2loop/thickness_calculator.py
@@ -782,7 +782,7 @@ def __init__(
             if sections is not None
             else geopandas.GeoDataFrame({"geometry": []}, geometry="geometry")  # ensure predictable structure when no sections supplied
         )
-        self.section_thickness_records: geopandas.GeoDataFrame | List[dict] = []  # populated as GeoDataFrame during compute
+        self.section_thickness_records: geopandas.GeoDataFrame = geopandas.GeoDataFrame()  # populated as GeoDataFrame during compute
         self.section_intersection_points: dict = {}
 
     def type(self):
@@ -894,7 +894,8 @@ def compute(
 
         try:
             geology_sindex = geology.sindex
-        except Exception:
+        except Exception as e:
+            logger.error("Failed to create spatial index for geology data: %s", e, exc_info=True)
             geology_sindex = None
 
         dip_column = next((col for col in ("DIP", "dip", "Dip") if col in structure_data.columns), None)
@@ -910,14 +911,19 @@ def compute(
                     try:
                         orientation_coords = orient_df[["X", "Y"]].astype(float).to_numpy()
                     except ValueError:
+                        logger.debug(
+                            "Failed to convert orientation coordinates to float for %d rows. Data quality issue likely. Example rows: %s",
+                            len(orient_df),
+                            orient_df[["X", "Y"]].head(3).to_dict(orient="records")
+                        )
                         orientation_coords = numpy.empty((0, 2))
                     if orientation_coords.size:
                         orientation_dips = orient_df["_dip_value"].astype(float).to_numpy()
                         try:
                             orientation_tree = cKDTree(orientation_coords)
-                        except Exception:
+                        except (ValueError, TypeError) as e:
+                            logger.error(f"Failed to construct cKDTree for orientation data: {e}")
                             orientation_tree = None
-
         default_dip_warning_emitted = False
 
         units_lookup = dict(zip(units["name"], units.index))
@@ -940,7 +946,8 @@ def compute(
                 continue
 
             split_segments = []
-            for _, poly in geology.iloc[candidate_idx].iterrows():
+            for idx in candidate_idx:
+                poly = geology.iloc[idx]
                 polygon_geom = poly.geometry
                 if polygon_geom is None or polygon_geom.is_empty:
                     continue
diff --git a/map2loop/utils.py b/map2loop/utils.py
index 927d4c5e..1434562b 100644
--- a/map2loop/utils.py
+++ b/map2loop/utils.py
@@ -629,7 +629,8 @@ def clean_line_geometry(geometry: shapely.geometry.base.BaseGeometry):
         return geometry
     try:
         merged = shapely.ops.linemerge(geometry)
-    except Exception:
+    except Exception as e:
+        logger.exception("Exception occurred in shapely.ops.linemerge during clean_line_geometry")
         return None
     if merged.geom_type == "LineString":
         return merged
@@ -737,7 +738,8 @@ def nearest_orientation_to_line(orientation_tree, orientation_dips, orientation_
     k = min(len(orientation_dips), 5)
     try:
         distances, indices = orientation_tree.query(query_xy, k=k)
-    except Exception:
+    except Exception as e:
+        logger.exception("Exception occurred during KDTree query in nearest_orientation_to_line; returning (nan, None, None)")
         return numpy.nan, None, None
     distances = numpy.atleast_1d(distances)
     indices = numpy.atleast_1d(indices)

From 6b8e67400b613e557afeb3c2de719f7deaa8515e Mon Sep 17 00:00:00 2001
From: lachlangrose <7371904+lachlangrose@users.noreply.github.com>
Date: Thu, 13 Nov 2025 01:35:34 +0000
Subject: [PATCH 05/31] style: style fixes by ruff and autoformatting by black

---
 map2loop/utils.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/map2loop/utils.py b/map2loop/utils.py
index 1434562b..3efe8e76 100644
--- a/map2loop/utils.py
+++ b/map2loop/utils.py
@@ -629,7 +629,7 @@ def clean_line_geometry(geometry: shapely.geometry.base.BaseGeometry):
         return geometry
     try:
         merged = shapely.ops.linemerge(geometry)
-    except Exception as e:
+    except Exception:
         logger.exception("Exception occurred in shapely.ops.linemerge during clean_line_geometry")
         return None
     if merged.geom_type == "LineString":
@@ -738,7 +738,7 @@ def nearest_orientation_to_line(orientation_tree, orientation_dips, orientation_
     k = min(len(orientation_dips), 5)
     try:
         distances, indices = orientation_tree.query(query_xy, k=k)
-    except Exception as e:
+    except Exception:
         logger.exception("Exception occurred during KDTree query in nearest_orientation_to_line; returning (nan, None, None)")
         return numpy.nan, None, None
     distances = numpy.atleast_1d(distances)

From 7a7ecbd22efa943a9e45983bc33f9f4e5f24596c Mon Sep 17 00:00:00 2001
From: rabii-chaarani <rabii.chaarani@monash.edu>
Date: Thu, 22 Jan 2026 10:19:01 +0930
Subject: [PATCH 06/31] refactor: added overlap checks

---
 map2loop/thickness_calculator.py | 59 +++++++++++++++++++++-----------
 1 file changed, 39 insertions(+), 20 deletions(-)

diff --git a/map2loop/thickness_calculator.py b/map2loop/thickness_calculator.py
index 614344cb..f6ca448c 100644
--- a/map2loop/thickness_calculator.py
+++ b/map2loop/thickness_calculator.py
@@ -862,25 +862,6 @@ def compute(
             )
             return units
 
-        unit_column_candidates = [
-            "UNITNAME",
-            "unitname",
-            "UNIT_NAME",
-            "UnitName",
-            "unit",
-            "Unit",
-            "UNIT",
-            "name",
-            "Name",
-        ]
-        unit_column = next((col for col in unit_column_candidates if col in geology_data.columns), None)
-        if unit_column is None:
-            logger.warning(
-                "AlongSection: Unable to identify a unit-name column in geology data; expected one of %s.",
-                unit_column_candidates,
-            )
-            return units
-
         sections = self.sections.copy()
         geology = geology_data.copy()
 
@@ -898,6 +879,44 @@ def compute(
             logger.error("Failed to create spatial index for geology data: %s", e, exc_info=True)
             geology_sindex = None
 
+        #TODO: check if sections and geology have same crs and reproject if needed. then check if sections overlap geology
+        sections_bounds = sections.total_bounds
+        geology_bounds = geology.total_bounds
+        if (
+            sections_bounds[2] < geology_bounds[0]
+            or sections_bounds[0] > geology_bounds[2]
+            or sections_bounds[3] < geology_bounds[1]
+            or sections_bounds[1] > geology_bounds[3]
+        ):
+            logger.warning("AlongSection: Sections do not overlap geology extent; skipping thickness calculation.")
+            return units
+
+        if geology_sindex is None:
+            overlaps = False
+            for geom in sections.geometry:
+                if geom is None or geom.is_empty:
+                    continue
+                if geology.intersects(geom).any():
+                    overlaps = True
+                    break
+            if not overlaps:
+                logger.warning("AlongSection: Sections do not overlap geology; skipping thickness calculation.")
+                return units
+        else:
+            overlaps = False
+            for geom in sections.geometry:
+                if geom is None or geom.is_empty:
+                    continue
+                candidate_idx = list(geology_sindex.intersection(geom.bounds))
+                if not candidate_idx:
+                    continue
+                if geology.iloc[candidate_idx].intersects(geom).any():
+                    overlaps = True
+                    break
+            if not overlaps:
+                logger.warning("AlongSection: Sections do not overlap geology; skipping thickness calculation.")
+                return units
+
         dip_column = next((col for col in ("DIP", "dip", "Dip") if col in structure_data.columns), None)
         orientation_tree = None
         orientation_dips = None
@@ -1006,7 +1025,7 @@ def compute(
                         )
                         default_dip_warning_emitted = True
 
-                thickness = segment["length"] * abs(math.sin(math.radians(dip_value)))
+                thickness = segment["length"] * math.sin(math.radians(dip_value))
                 if thickness <= 0:
                     continue
 

From d8bc0609d6f6fa93c609f945e7757e642fc2f565 Mon Sep 17 00:00:00 2001
From: rabii-chaarani <rabii.chaarani@monash.edu>
Date: Thu, 22 Jan 2026 10:22:45 +0930
Subject: [PATCH 07/31] fix: use standard unit column name

---
 map2loop/thickness_calculator.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/map2loop/thickness_calculator.py b/map2loop/thickness_calculator.py
index f6ca448c..83f2a290 100644
--- a/map2loop/thickness_calculator.py
+++ b/map2loop/thickness_calculator.py
@@ -861,7 +861,7 @@ def compute(
                 "AlongSection: Geology polygons are required to split sections; skipping thickness calculation."
             )
             return units
-
+        unit_column = "UNITNAME"
         sections = self.sections.copy()
         geology = geology_data.copy()
 

From 431bc0a96bbdd10084c645c8df04cd8caa7fd7d6 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Mon, 1 Dec 2025 16:04:23 +1100
Subject: [PATCH 08/31] fix: move extra sort arguments to sorter init (#241)

* fix: move extra arguments to init

* Apply suggestions from code review

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>

* fix: resolving copilot review and adding class attribute with required arguments

* fix: remove SorterUseHint call

* fix: remove arguments from ABC
make all requrements a list
make all arguments required

* style: style fixes by ruff and autoformatting by black

* fix: add check for geology_data

* fix: not initiliase sorter in project.init

* fix: try use SorterUseNetworkX

* fix: adding utils module and fixing tests

* style: style fixes by ruff and autoformatting by black

* fix: sorter arguments optional. Allow project to update attributes from map data

* fix: minage/maxage optional

* fix: add unit relations as optional argument

* fix: default min/max age

* fix: wrong defaults.

---------

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
Co-authored-by: lachlangrose <7371904+lachlangrose@users.noreply.github.com>
Co-authored-by: rabii-chaarani <rabii.chaarani@monash.edu>
---
 map2loop/__init__.py                          |   6 +-
 map2loop/logging.py                           |  18 +-
 map2loop/project.py                           |  33 +-
 map2loop/sorter.py                            | 289 ++++++++++--------
 map2loop/utils/__init__.py                    |  18 ++
 map2loop/utils/load_map2loop_data.py          |  21 ++
 .../{utils.py => utils/utility_functions.py}  |   2 +-
 tests/project/test_plot_hamersley.py          |   2 +
 tests/project/test_thickness_calculations.py  |   4 +-
 .../test_interpolated_structure.py            |   2 +-
 .../test_ThicknessStructuralPoint.py          |   2 +-
 .../test_ThicknessCalculatorAlpha.py          |   2 +-
 12 files changed, 244 insertions(+), 155 deletions(-)
 create mode 100644 map2loop/utils/__init__.py
 create mode 100644 map2loop/utils/load_map2loop_data.py
 rename map2loop/{utils.py => utils/utility_functions.py} (99%)

diff --git a/map2loop/__init__.py b/map2loop/__init__.py
index 8723f4ef..66a14c9a 100644
--- a/map2loop/__init__.py
+++ b/map2loop/__init__.py
@@ -1,10 +1,6 @@
 import logging
+from map2loop.logging import loggers, ch
 
-loggers = {}
-ch = logging.StreamHandler()
-formatter = logging.Formatter("%(levelname)s: %(asctime)s: %(filename)s:%(lineno)d -- %(message)s")
-ch.setFormatter(formatter)
-ch.setLevel(logging.WARNING)
 from .project import Project
 from .version import __version__
 
diff --git a/map2loop/logging.py b/map2loop/logging.py
index 2daaa0c2..816ac45e 100644
--- a/map2loop/logging.py
+++ b/map2loop/logging.py
@@ -1,7 +1,11 @@
 import logging
-import map2loop
 
+loggers = {}
 
+ch = ch = logging.StreamHandler()
+formatter = logging.Formatter("%(levelname)s: %(asctime)s: %(filename)s:%(lineno)d -- %(message)s")
+ch.setFormatter(formatter)
+ch.setLevel(logging.WARNING)
 def get_levels():
     """dict for converting to logger levels from string
 
@@ -33,12 +37,12 @@ def getLogger(name: str):
     logging.Logger
         logger object
     """
-    if name in map2loop.loggers:
-        return map2loop.loggers[name]
+    if name in loggers:
+        return loggers[name]
     logger = logging.getLogger(name)
-    logger.addHandler(map2loop.ch)
+    logger.addHandler(ch)
     logger.propagate = False
-    map2loop.loggers[name] = logger
+    loggers[name] = logger
     return logger
 
 
@@ -56,9 +60,9 @@ def set_level(level: str):
     """
     levels = get_levels()
     level = levels.get(level, logging.WARNING)
-    map2loop.ch.setLevel(level)
+    ch.setLevel(level)
 
-    for name in map2loop.loggers:
+    for name in loggers:
         logger = logging.getLogger(name)
         logger.setLevel(level)
     logger.info(f"Logging level set to {level}")
diff --git a/map2loop/project.py b/map2loop/project.py
index cf5840d1..ec7f10f5 100644
--- a/map2loop/project.py
+++ b/map2loop/project.py
@@ -140,7 +140,7 @@ def __init__(
         self.set_default_samplers()
         self.bounding_box = bounding_box
         self.contact_extractor = None
-        self.sorter = SorterUseHint()
+        self.sorter = SorterUseHint
         self.throw_calculator = ThrowCalculatorAlpha()
         self.fault_orientation = FaultOrientationNearest()
         self.map_data = MapData(verbose_level=verbose_level)
@@ -560,7 +560,15 @@ def calculate_stratigraphic_order(self, take_best=False):
             )
             self.contact_extractor.extract_all_contacts()
         if take_best:
-            sorters = [SorterUseHint(), SorterAgeBased(), SorterAlpha(), SorterUseNetworkX()]
+            sorters = [
+                SorterAgeBased(),
+                SorterAlpha(
+                    contacts=self.contact_extractor.contacts,
+                ),
+                SorterUseNetworkX(
+                    unit_relationships=self.topology.get_unit_unit_relationships(),
+                ),
+            ]
             logger.info(
                 f"Calculating best stratigraphic column from {[sorter.sorter_label for sorter in sorters]}"
             )
@@ -568,11 +576,6 @@ def calculate_stratigraphic_order(self, take_best=False):
             columns = [
                 sorter.sort(
                     self.stratigraphic_column.stratigraphicUnits,
-                    self.topology.get_unit_unit_relationships(),
-                    self.contact_extractor.contacts,
-                    self.map_data.get_map_data(Datatype.GEOLOGY),
-                    self.map_data.get_map_data(Datatype.STRUCTURE),
-                    self.map_data.get_map_data(Datatype.DTM),
                 )
                 for sorter in sorters
             ]
@@ -600,13 +603,19 @@ def calculate_stratigraphic_order(self, take_best=False):
             self.stratigraphic_column.column = column
         else:
             logger.info(f'Calculating stratigraphic column using sorter {self.sorter.sorter_label}')
+            # Update sorter with current data based on what it needs
+            if hasattr(self.sorter, 'unit_relationships') and self.sorter.unit_relationships is None:
+                self.sorter.unit_relationships = self.topology.get_unit_unit_relationships()
+            if hasattr(self.sorter, 'contacts') and self.sorter.contacts is None:
+                self.sorter.contacts = self.contact_extractor.contacts
+            if hasattr(self.sorter, 'geology_data') and self.sorter.geology_data is None:
+                self.sorter.geology_data = self.map_data.get_map_data(Datatype.GEOLOGY)
+            if hasattr(self.sorter, 'structure_data') and self.sorter.structure_data is None:
+                self.sorter.structure_data = self.map_data.get_map_data(Datatype.STRUCTURE)
+            if hasattr(self.sorter, 'dtm_data') and self.sorter.dtm_data is None:
+                self.sorter.dtm_data = self.map_data.get_map_data(Datatype.DTM)
             self.stratigraphic_column.column = self.sorter.sort(
                 self.stratigraphic_column.stratigraphicUnits,
-                self.topology.get_unit_unit_relationships(),
-                self.contact_extractor.contacts,
-                self.map_data.get_map_data(Datatype.GEOLOGY),
-                self.map_data.get_map_data(Datatype.STRUCTURE),
-                self.map_data.get_map_data(Datatype.DTM),
             )
 
     @beartype.beartype
diff --git a/map2loop/sorter.py b/map2loop/sorter.py
index 0428ace6..7d1a7672 100644
--- a/map2loop/sorter.py
+++ b/map2loop/sorter.py
@@ -3,10 +3,10 @@
 import pandas
 import numpy as np
 import math
-from typing import Union
-from osgeo import gdal
+from typing import Union, Optional, List
+from map2loop.topology import Topology
 import geopandas
-
+from osgeo import gdal
 from map2loop.utils import value_from_raster
 from .logging import getLogger
 
@@ -21,11 +21,20 @@ class Sorter(ABC):
         ABC (ABC): Derived from Abstract Base Class
     """
 
-    def __init__(self):
+    def __init__(
+        self):
         """
-        Initialiser of for Sorter
+        Initialiser for Sorter
+        
+        Args:
+            unit_relationships (pandas.DataFrame): the relationships between units (columns must contain ["Index1", "Unitname1", "Index2", "Unitname2"])
+            contacts (pandas.DataFrame): unit contacts with length of the contacts in metres
+            geology_data (geopandas.GeoDataFrame): the geology data
+            structure_data (geopandas.GeoDataFrame): the structure data
+            dtm_data (gdal.Dataset): the dtm data
         """
         self.sorter_label = "SorterBaseClass"
+        
 
     def type(self):
         """
@@ -38,25 +47,12 @@ def type(self):
 
     @beartype.beartype
     @abstractmethod
-    def sort(
-        self,
-        units: pandas.DataFrame,
-        unit_relationships: pandas.DataFrame,
-        contacts: pandas.DataFrame,
-        geology_data: geopandas.GeoDataFrame = None,
-        structure_data: geopandas.GeoDataFrame = None,
-        dtm_data: gdal.Dataset = None,
-    ) -> list:
+    def sort(self, units: pandas.DataFrame) -> list:
         """
         Execute sorter method (abstract method)
 
         Args:
             units (pandas.DataFrame): the data frame to sort (columns must contain ["layerId", "name", "minAge", "maxAge", "group"])
-            units_relationships (pandas.DataFrame): the relationships between units (columns must contain ["Index1", "Unitname1", "Index2", "Unitname2"])
-            contacts (pandas.DataFrame): unit contacts with length of the contacts in metres
-            geology_data (geopandas.GeoDataFrame): the geology data
-            structure_data (geopandas.GeoDataFrame): the structure data
-            dtm_data (ggdal.Dataset): the dtm data
 
         Returns:
             list: sorted list of unit names
@@ -68,42 +64,71 @@ class SorterUseNetworkX(Sorter):
     """
     Sorter class which returns a sorted list of units based on the unit relationships using a topological graph sorting algorithm
     """
+    required_arguments: List[str] = [
+        'geology_data'
+    ]
 
-    def __init__(self):
+    def __init__(
+        self,
+        *,
+        unit_relationships: Optional[pandas.DataFrame] = None,
+        geology_data: Optional[geopandas.GeoDataFrame] = None,
+    ):
         """
         Initialiser for networkx graph sorter
+        
+        Args:
+            unit_relationships (pandas.DataFrame): the relationships between units
         """
+        super().__init__()
         self.sorter_label = "SorterUseNetworkX"
+        if geology_data is not None:
+            self.set_geology_data(geology_data)
+        elif unit_relationships is not None:
+            self.unit_relationships = unit_relationships
+        else:
+            self.unit_relationships = None
+    def set_geology_data(self, geology_data: geopandas.GeoDataFrame):
+        """
+        Set geology data and calculate topology and unit relationships
 
+        Args:
+            geology_data (geopandas.GeoDataFrame): the geology data
+        """
+        self._calculate_topology(geology_data)
+    def _calculate_topology(self, geology_data: geopandas.GeoDataFrame):
+        if not geology_data:
+            raise ValueError("geology_data is required")
+
+        if isinstance(geology_data, geopandas.GeoDataFrame) is False:
+            raise TypeError("geology_data must be a geopandas.GeoDataFrame")
+        
+        if 'UNITNAME' not in geology_data.columns:
+            raise ValueError("geology_data must contain 'UNITNAME' column")
+        
+        self.topology = Topology(geology_data=geology_data)
+        self.unit_relationships = self.topology.get_unit_unit_relationships()
+    
     @beartype.beartype
-    def sort(
-        self,
-        units: pandas.DataFrame,
-        unit_relationships: pandas.DataFrame,
-        contacts: pandas.DataFrame,
-        geology_data: geopandas.GeoDataFrame = None,
-        structure_data: geopandas.GeoDataFrame = None,
-        dtm_data: gdal.Dataset = None,
-    ) -> list:
+    def sort(self, units: pandas.DataFrame) -> list:
         """
-        Execute sorter method takes unit data, relationships and a hint and returns the sorted unit names based on this algorithm.
+        Execute sorter method takes unit data and returns the sorted unit names based on this algorithm.
 
         Args:
             units (pandas.DataFrame): the data frame to sort
-            units_relationships (pandas.DataFrame): the relationships between units
-            contacts (pandas.DataFrame): unit contacts with length of the contacts in metres
 
         Returns:
             list: the sorted unit names
         """
         import networkx as nx
-
+        if self.unit_relationships is None:
+            raise ValueError("SorterUseNetworkX requires 'unit_relationships' argument")
         graph = nx.DiGraph()
         name_to_index = {}
         for row in units.iterrows():
             graph.add_node(int(row[1]["layerId"]), name=row[1]["name"])
             name_to_index[row[1]["name"]] = int(row[1]["layerId"])
-        for row in unit_relationships.iterrows():
+        for row in self.unit_relationships.iterrows():
             graph.add_edge(name_to_index[row[1]["UNITNAME_1"]], name_to_index[row[1]["UNITNAME_2"]])
 
         cycles = list(nx.simple_cycles(graph))
@@ -124,54 +149,55 @@ def sort(
 
 
 class SorterUseHint(SorterUseNetworkX):
-    def __init__(self):
-        logger.info(
-            "SorterUseHint is deprecated in v3.2. Use SorterUseNetworkX instead"
+    required_arguments: List[str] = ['unit_relationships']
+    def __init__(
+        self,
+        *,
+        geology_data: Optional[geopandas.GeoDataFrame] = None,
+    ):
+        logger.warning(
+            "SorterUseHint is deprecated in v3.2. Using SorterUseNetworkX instead"
         )
-        super().__init__()
+        super().__init__(geology_data=geology_data)
+    
 
 
 class SorterAgeBased(Sorter):
     """
     Sorter class which returns a sorted list of units based on the min and max ages of the units
     """
-
-    def __init__(self):
+    required_arguments = ['min_age_column','max_age_column']
+    def __init__(self,*, min_age_column:Optional[str]='minAge', max_age_column:Optional[str]='maxAge'):
         """
         Initialiser for age based sorter
         """
+        super().__init__()
+        self.min_age_column = min_age_column
+        self.max_age_column = max_age_column
         self.sorter_label = "SorterAgeBased"
 
-    def sort(
-        self,
-        units: pandas.DataFrame,
-        unit_relationships: pandas.DataFrame,
-        contacts: pandas.DataFrame,
-        geology_data: geopandas.GeoDataFrame = None,
-        structure_data: geopandas.GeoDataFrame = None,
-        dtm_data: gdal.Dataset = None,
-    ) -> list:
+    def sort(self, units: pandas.DataFrame) -> list:
         """
-        Execute sorter method takes unit data, relationships and a hint and returns the sorted unit names based on this algorithm.
+        Execute sorter method takes unit data and returns the sorted unit names based on this algorithm.
 
         Args:
             units (pandas.DataFrame): the data frame to sort
-            units_relationships (pandas.DataFrame): the relationships between units
-            stratigraphic_order_hint (list): a list of unit names to use as a hint to sorting the units
-            contacts (pandas.DataFrame): unit contacts with length of the contacts in metres
 
         Returns:
             list: the sorted unit names
         """
         logger.info("Calling age based sorter")
         sorted_units = units.copy()
-        if "minAge" in units.columns and "maxAge" in units.columns:
+
+        if self.min_age_column in units.columns and self.max_age_column in units.columns:
             # print(sorted_units["minAge"], sorted_units["maxAge"])
             sorted_units["meanAge"] = sorted_units.apply(
-                lambda row: (row["minAge"] + row["maxAge"]) / 2.0, axis=1
+                lambda row: (row[self.min_age_column] + row[self.max_age_column]) / 2.0, axis=1
             )
         else:
-            sorted_units["meanAge"] = 0
+            logger.error(f"Columns {self.min_age_column} and {self.max_age_column} must be present in units DataFrame")
+            logger.error(f"Available columns are: {units.columns.tolist()}")
+            raise ValueError(f"Columns {self.min_age_column} and {self.max_age_column} must be present in units DataFrame")
         if "group" in units.columns:
             sorted_units = sorted_units.sort_values(by=["group", "meanAge"])
         else:
@@ -188,45 +214,48 @@ class SorterAlpha(Sorter):
     Sorter class which returns a sorted list of units based on the adjacency of units
     prioritising the units with lower number of contacting units
     """
-
-    def __init__(self):
+    required_arguments = ['contacts']
+    def __init__(
+        self,
+        *,
+        contacts: Optional[geopandas.GeoDataFrame] = None,
+    ):
         """
         Initialiser for adjacency based sorter
+        
+        Args:
+            contacts (geopandas.GeoDataFrame): unit contacts with length of the contacts in metres
         """
+        super().__init__()
+        self.contacts = contacts
         self.sorter_label = "SorterAlpha"
+        if 'UNITNAME_1' not in contacts.columns or 'UNITNAME_2' not in contacts.columns or 'length' not in contacts.columns:
+            raise ValueError("contacts GeoDataFrame must contain 'UNITNAME_1', 'UNITNAME_2' and 'length' columns")
 
-    def sort(
-        self,
-        units: pandas.DataFrame,
-        unit_relationships: pandas.DataFrame,
-        contacts: pandas.DataFrame,
-        geology_data: geopandas.GeoDataFrame = None,
-        structure_data: geopandas.GeoDataFrame = None,
-        dtm_data: gdal.Dataset = None,
-    ) -> list:
+    def sort(self, units: pandas.DataFrame) -> list:
         """
-        Execute sorter method takes unit data, relationships and a hint and returns the sorted unit names based on this algorithm.
+        Execute sorter method takes unit data and returns the sorted unit names based on this algorithm.
 
         Args:
             units (pandas.DataFrame): the data frame to sort
-            units_relationships (pandas.DataFrame): the relationships between units
-            stratigraphic_order_hint (list): a list of unit names to use as a hint to sorting the units
-            contacts (pandas.DataFrame): unit contacts with length of the contacts in metres
 
         Returns:
             list: the sorted unit names
         """
+        if self.contacts is None:
+            raise ValueError("contacts must be set (not None) before calling sort() in SorterAlpha.")
         import networkx as nx
-
-        contacts = contacts.sort_values(by="length", ascending=False)[
+        if self.contacts is None:
+            raise ValueError("SorterAlpha requires 'contacts' argument")
+        sorted_contacts = self.contacts.sort_values(by="length", ascending=False)[
             ["UNITNAME_1", "UNITNAME_2", "length"]
         ]
-        units = list(units["name"].unique())
+        unit_names = list(units["name"].unique())
         graph = nx.Graph()
-        for unit in units:
+        for unit in unit_names:
             graph.add_node(unit, name=unit)
-        max_weight = max(list(contacts["length"])) + 1
-        for _, row in contacts.iterrows():
+        max_weight = max(list(sorted_contacts["length"])) + 1
+        for _, row in sorted_contacts.iterrows():
             graph.add_edge(
                 row["UNITNAME_1"], row["UNITNAME_2"], weight=int(max_weight - row["length"])
             )
@@ -272,51 +301,51 @@ class SorterMaximiseContacts(Sorter):
     Sorter class which returns a sorted list of units based on the adjacency of units
     prioritising the maximum length of each contact
     """
-
-    def __init__(self):
+    required_arguments = ['contacts']
+    def __init__(
+        self,
+        *,
+        contacts: Optional[geopandas.GeoDataFrame] = None,
+    ):
         """
         Initialiser for adjacency based sorter
-
+        
+        Args:
+            contacts (pandas.DataFrame): unit contacts with length of the contacts in metres
         """
+        super().__init__()
         self.sorter_label = "SorterMaximiseContacts"
         # variables for visualising/interrogating the sorter
         self.graph = None
         self.route = None
         self.directed_graph = None
+        self.contacts = contacts
+        if 'UNITNAME_1' not in contacts.columns or 'UNITNAME_2' not in contacts.columns or 'length' not in contacts.columns:
+            raise ValueError("contacts GeoDataFrame must contain 'UNITNAME_1', 'UNITNAME_2' and 'length' columns")
 
-    def sort(
-        self,
-        units: pandas.DataFrame,
-        unit_relationships: pandas.DataFrame,
-        contacts: pandas.DataFrame,
-        geology_data: geopandas.GeoDataFrame = None,
-        structure_data: geopandas.GeoDataFrame = None,
-        dtm_data: gdal.Dataset = None,
-    ) -> list:
+    def sort(self, units: pandas.DataFrame) -> list:
         """
-        Execute sorter method takes unit data, relationships and a hint and returns the sorted unit names based on this algorithm.
+        Execute sorter method takes unit data and returns the sorted unit names based on this algorithm.
 
         Args:
             units (pandas.DataFrame): the data frame to sort
-            units_relationships (pandas.DataFrame): the relationships between units
-            stratigraphic_order_hint (list): a list of unit names to use as a hint to sorting the units
-            contacts (pandas.DataFrame): unit contacts with length of the contacts in metres
 
         Returns:
             list: the sorted unit names
         """
         import networkx as nx
         import networkx.algorithms.approximation as nx_app
-
-        sorted_contacts = contacts.sort_values(by="length", ascending=False)
+        if self.contacts is None:
+            raise ValueError("SorterMaximiseContacts requires 'contacts' argument")
+        sorted_contacts = self.contacts.sort_values(by="length", ascending=False)
         self.graph = nx.Graph()
-        units = list(units["name"].unique())
-        for unit in units:
+        unit_names = list(units["name"].unique())
+        for unit in unit_names:
             ## some units may not have any contacts e.g. if they are intrusives or sills. If we leave this then the
             ## sorter crashes
             if (
-                unit not in sorted_contacts['UNITNAME_1']
-                or unit not in sorted_contacts['UNITNAME_2']
+                unit not in sorted_contacts['UNITNAME_1'].values
+                or unit not in sorted_contacts['UNITNAME_2'].values
             ):
                 continue
             self.graph.add_node(unit, name=unit)
@@ -355,38 +384,41 @@ class SorterObservationProjections(Sorter):
     Sorter class which returns a sorted list of units based on the adjacency of units
     using the direction of observations to predict which unit is adjacent to the current one
     """
-
-    def __init__(self, length: Union[float, int] = 1000):
+    required_arguments = ['contacts', 'geology_data', 'structure_data', 'dtm_data']
+    def __init__(
+        self,
+        *,
+        contacts: Optional[geopandas.GeoDataFrame] = None,
+        geology_data: Optional[geopandas.GeoDataFrame] = None,
+        structure_data: Optional[geopandas.GeoDataFrame] = None,
+        dtm_data: Optional[gdal.Dataset] = None,
+        length: Union[float, int] = 1000
+    ):
         """
         Initialiser for adjacency based sorter
 
         Args:
+            contacts (pandas.DataFrame): unit contacts with length of the contacts in metres
+            geology_data (geopandas.GeoDataFrame): the geology data
+            structure_data (geopandas.GeoDataFrame): the structure data
+            dtm_data (gdal.Dataset): the dtm data
             length (int): the length of the projection in metres
         """
+        super().__init__()
+        self.contacts = contacts
+        self.geology_data = geology_data
+        self.structure_data = structure_data
+        self.dtm_data = dtm_data
         self.sorter_label = "SorterObservationProjections"
         self.length = length
         self.lines = []
 
-    def sort(
-        self,
-        units: pandas.DataFrame,
-        unit_relationships: pandas.DataFrame,
-        contacts: pandas.DataFrame,
-        geology_data: geopandas.GeoDataFrame,
-        structure_data: geopandas.GeoDataFrame,
-        dtm_data: gdal.Dataset
-    ) -> list:
+    def sort(self, units: pandas.DataFrame) -> list:
         """
-        Execute sorter method takes unit data, relationships and a hint and returns the sorted unit names based on this algorithm.
+        Execute sorter method takes unit data and returns the sorted unit names based on this algorithm.
 
         Args:
             units (pandas.DataFrame): the data frame to sort
-            units_relationships (pandas.DataFrame): the relationships between units
-            stratigraphic_order_hint (list): a list of unit names to use as a hint to sorting the units
-            contacts (pandas.DataFrame): unit contacts with length of the contacts in metres
-            geology_data (geopandas.GeoDataFrame): the geology data
-            structure_data (geopandas.GeoDataFrame): the structure data
-            dtm_data (ggdal.Dataset): the dtm data
 
         Returns:
             list: the sorted unit names
@@ -394,15 +426,22 @@ def sort(
         import networkx as nx
         import networkx.algorithms.approximation as nx_app
         from shapely.geometry import LineString, Point
-
-        geol = geology_data.copy()
+        if self.contacts is None:
+            raise ValueError("SorterObservationProjections requires 'contacts' argument")
+        if self.geology_data is None:
+            raise ValueError("SorterObservationProjections requires 'geology_data' argument")
+        geol = self.geology_data.copy()
         if "INTRUSIVE" in geol.columns:
             geol = geol.drop(geol.index[geol["INTRUSIVE"]])
         if "SILL" in geol.columns:
             geol = geol.drop(geol.index[geol["SILL"]])
-        orientations = structure_data.copy()
-        inv_geotransform = gdal.InvGeoTransform(dtm_data.GetGeoTransform())
-        dtm_array = np.array(dtm_data.GetRasterBand(1).ReadAsArray().T)
+        if self.structure_data is None:
+            raise ValueError("structure_data is required for sorting but is None.")
+        orientations = self.structure_data.copy()
+        if self.dtm_data is None:
+            raise ValueError("DTM data (self.dtm_data) is not set. Cannot proceed with sorting.")
+        inv_geotransform = gdal.InvGeoTransform(self.dtm_data.GetGeoTransform())
+        dtm_array = np.array(self.dtm_data.GetRasterBand(1).ReadAsArray().T)
 
         # Create a map of maps to store younger/older observations
         ordered_unit_observations = []
@@ -504,9 +543,9 @@ def sort(
         g_undirected = g.to_undirected()
         for unit in unit_names:
             if len(list(g_undirected.neighbors(unit))) < 1:
-                mask1 = contacts["UNITNAME_1"] == unit
-                mask2 = contacts["UNITNAME_2"] == unit
-                for _, row in contacts[mask1 | mask2].iterrows():
+                mask1 = self.contacts["UNITNAME_1"] == unit
+                mask2 = self.contacts["UNITNAME_2"] == unit
+                for _, row in self.contacts[mask1 | mask2].iterrows():
                     if unit == row["UNITNAME_1"]:
                         g.add_edge(row["UNITNAME_2"], unit, weight=max_value * 10)
                     else:
diff --git a/map2loop/utils/__init__.py b/map2loop/utils/__init__.py
new file mode 100644
index 00000000..c7fa49d6
--- /dev/null
+++ b/map2loop/utils/__init__.py
@@ -0,0 +1,18 @@
+from .utility_functions import (
+    set_z_values_from_raster_df,
+    value_from_raster,
+    update_from_legacy_file,
+    preprocess_hjson_to_json,
+    read_hjson_with_json,
+    calculate_endpoints,
+    calculate_minimum_fault_length,
+    hex_to_rgb,
+    generate_random_hex_colors,
+    rebuild_sampled_basal_contacts,
+    multiline_to_line,
+    find_segment_strike_from_pt,
+    create_points,
+    normal_vector_to_dipdirection_dip,
+    strike_dip_vector,
+    generate_grid,
+)
diff --git a/map2loop/utils/load_map2loop_data.py b/map2loop/utils/load_map2loop_data.py
new file mode 100644
index 00000000..6cedb18b
--- /dev/null
+++ b/map2loop/utils/load_map2loop_data.py
@@ -0,0 +1,21 @@
+import geopandas
+import pathlib
+from osgeo import gdal
+gdal.UseExceptions()
+
+# Use the path of this file to locate the datasets directory
+def map2loop_dir(folder) -> pathlib.Path:
+    path = pathlib.Path(__file__).parent.parent.parent / 'map2loop' / '_datasets' / 'geodata_files' / f'{folder}'
+    return path
+
+def load_hamersley_geology():
+    path = map2loop_dir('hamersley') / "geology.geojson"
+    return geopandas.read_file(str(path))
+
+def load_hamersley_structure():
+    path = map2loop_dir('hamersley') / "structure.geojson"
+    return geopandas.read_file(str(path))
+
+def load_hamersley_dtm():
+    path = map2loop_dir('hamersley') / "dtm_rp.tif"
+    return gdal.Open(str(path))
diff --git a/map2loop/utils.py b/map2loop/utils/utility_functions.py
similarity index 99%
rename from map2loop/utils.py
rename to map2loop/utils/utility_functions.py
index 3efe8e76..082c615d 100644
--- a/map2loop/utils.py
+++ b/map2loop/utils/utility_functions.py
@@ -9,7 +9,7 @@
 import json
 from osgeo import gdal
 
-from .logging import getLogger
+from ..logging import getLogger
 logger = getLogger(__name__)
 
 
diff --git a/tests/project/test_plot_hamersley.py b/tests/project/test_plot_hamersley.py
index 504c4585..dadbb7d7 100644
--- a/tests/project/test_plot_hamersley.py
+++ b/tests/project/test_plot_hamersley.py
@@ -1,5 +1,6 @@
 import pytest
 from map2loop.project import Project
+from map2loop.sorter import SorterUseNetworkX
 from map2loop.m2l_enums import VerboseLevel
 from unittest.mock import patch
 from pyproj.exceptions import CRSError
@@ -36,6 +37,7 @@ def test_project_execution():
     except Exception:
         pytest.skip("Skipping the project test from server data due to loading failure")
     try:
+        proj.set_sorter(SorterUseNetworkX())
         proj.run_all(take_best=True)
     except requests.exceptions.ReadTimeout:
         pytest.skip(
diff --git a/tests/project/test_thickness_calculations.py b/tests/project/test_thickness_calculations.py
index 373687ed..961e1a9c 100644
--- a/tests/project/test_thickness_calculations.py
+++ b/tests/project/test_thickness_calculations.py
@@ -3,11 +3,10 @@
 import geopandas
 import numpy
 
-from map2loop._datasets.geodata_files import load_map2loop_data
+from map2loop.utils import load_map2loop_data
 from map2loop.thickness_calculator import InterpolatedStructure, StructuralPoint
 from map2loop import Project
 
-
 # 1. self.stratigraphic_column.stratigraphicUnits,
 st_units = pandas.DataFrame(
     {
@@ -1705,6 +1704,7 @@ def test_calculate_unit_thicknesses():
     ], "Default for thickness calculator not set"  ## default is InterpolatedStructure
 
     # check set
+    print("****",StructuralPoint.__module__)
 
     project.set_thickness_calculator([StructuralPoint(dtm_data=dtm, bounding_box=bbox_3d), InterpolatedStructure(dtm_data=dtm, bounding_box=bbox_3d)])
     assert project.get_thickness_calculator() == [
diff --git a/tests/thickness/InterpolatedStructure/test_interpolated_structure.py b/tests/thickness/InterpolatedStructure/test_interpolated_structure.py
index 164ce2f7..bbda4e35 100644
--- a/tests/thickness/InterpolatedStructure/test_interpolated_structure.py
+++ b/tests/thickness/InterpolatedStructure/test_interpolated_structure.py
@@ -4,7 +4,7 @@
 
 from map2loop.mapdata import MapData
 from map2loop.thickness_calculator import InterpolatedStructure
-from map2loop._datasets.geodata_files.load_map2loop_data import (
+from map2loop.utils.load_map2loop_data import (
     load_hamersley_geology,
     load_hamersley_dtm,
 )
diff --git a/tests/thickness/StructurePoint/test_ThicknessStructuralPoint.py b/tests/thickness/StructurePoint/test_ThicknessStructuralPoint.py
index 74c7851d..0fbc1185 100644
--- a/tests/thickness/StructurePoint/test_ThicknessStructuralPoint.py
+++ b/tests/thickness/StructurePoint/test_ThicknessStructuralPoint.py
@@ -4,7 +4,7 @@
 
 from map2loop.mapdata import MapData
 from map2loop.thickness_calculator import StructuralPoint
-from map2loop._datasets.geodata_files.load_map2loop_data import load_hamersley_geology
+from map2loop.utils.load_map2loop_data import load_hamersley_geology
 from map2loop.m2l_enums import Datatype
 
 ####################################################################
diff --git a/tests/thickness/ThicknessCalculatorAlpha/test_ThicknessCalculatorAlpha.py b/tests/thickness/ThicknessCalculatorAlpha/test_ThicknessCalculatorAlpha.py
index 521b633f..ba740a18 100644
--- a/tests/thickness/ThicknessCalculatorAlpha/test_ThicknessCalculatorAlpha.py
+++ b/tests/thickness/ThicknessCalculatorAlpha/test_ThicknessCalculatorAlpha.py
@@ -5,7 +5,7 @@
 from map2loop.mapdata import MapData
 from map2loop.m2l_enums import Datatype
 from map2loop.thickness_calculator import ThicknessCalculatorAlpha
-from map2loop._datasets.geodata_files.load_map2loop_data import load_hamersley_geology
+from map2loop.utils.load_map2loop_data import load_hamersley_geology
 
 
 #########################################################

From 1451d36d427c7526d6faa14a3e8f1275489c552c Mon Sep 17 00:00:00 2001
From: "github-actions[bot]"
 <41898282+github-actions[bot]@users.noreply.github.com>
Date: Mon, 1 Dec 2025 15:00:56 +0930
Subject: [PATCH 09/31] chore(master): release 3.2.9 (#242)

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
---
 .release-please-manifest.json | 2 +-
 CHANGELOG.md                  | 7 +++++++
 map2loop/version.py           | 2 +-
 3 files changed, 9 insertions(+), 2 deletions(-)

diff --git a/.release-please-manifest.json b/.release-please-manifest.json
index c4048623..bbadce2a 100644
--- a/.release-please-manifest.json
+++ b/.release-please-manifest.json
@@ -1,3 +1,3 @@
 {
-  ".": "3.2.8"
+  ".": "3.2.9"
 }
\ No newline at end of file
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9ed7373d..2729cc98 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,12 @@
 # Changelog
 
+## [3.2.9](https://github.com/Loop3D/map2loop/compare/v3.2.8...v3.2.9) (2025-12-01)
+
+
+### Bug Fixes
+
+* move extra sort arguments to sorter init ([#241](https://github.com/Loop3D/map2loop/issues/241)) ([3414fec](https://github.com/Loop3D/map2loop/commit/3414fec29969407f7b1b07ee59bfbedfe90a63b1))
+
 ## [3.2.8](https://github.com/Loop3D/map2loop/compare/v3.2.7...v3.2.8) (2025-10-29)
 
 
diff --git a/map2loop/version.py b/map2loop/version.py
index 78e82ed1..b57bdcdb 100644
--- a/map2loop/version.py
+++ b/map2loop/version.py
@@ -1 +1 @@
-__version__ = "3.2.8"
+__version__ = "3.2.9"

From c60264de6c6b815ba16fde4f7c6264440674b45f Mon Sep 17 00:00:00 2001
From: lachlangrose <lachlan.grose@monash.edu>
Date: Tue, 2 Dec 2025 10:48:56 +1100
Subject: [PATCH 10/31] fix: change to no default arguments and add
 unit_name_column as a required

---
 map2loop/sorter.py | 58 ++++++++++++++++++++++++++++++----------------
 1 file changed, 38 insertions(+), 20 deletions(-)

diff --git a/map2loop/sorter.py b/map2loop/sorter.py
index 7d1a7672..738cfd25 100644
--- a/map2loop/sorter.py
+++ b/map2loop/sorter.py
@@ -65,12 +65,14 @@ class SorterUseNetworkX(Sorter):
     Sorter class which returns a sorted list of units based on the unit relationships using a topological graph sorting algorithm
     """
     required_arguments: List[str] = [
-        'geology_data'
+        'geology_data',
+        'unit_name_column'
     ]
 
     def __init__(
         self,
         *,
+        unit_name_column:Optional[str]='name',
         unit_relationships: Optional[pandas.DataFrame] = None,
         geology_data: Optional[geopandas.GeoDataFrame] = None,
     ):
@@ -82,6 +84,7 @@ def __init__(
         """
         super().__init__()
         self.sorter_label = "SorterUseNetworkX"
+        self.unit_name_column = unit_name_column
         if geology_data is not None:
             self.set_geology_data(geology_data)
         elif unit_relationships is not None:
@@ -97,7 +100,7 @@ def set_geology_data(self, geology_data: geopandas.GeoDataFrame):
         """
         self._calculate_topology(geology_data)
     def _calculate_topology(self, geology_data: geopandas.GeoDataFrame):
-        if not geology_data:
+        if geology_data is None:
             raise ValueError("geology_data is required")
 
         if isinstance(geology_data, geopandas.GeoDataFrame) is False:
@@ -166,12 +169,13 @@ class SorterAgeBased(Sorter):
     """
     Sorter class which returns a sorted list of units based on the min and max ages of the units
     """
-    required_arguments = ['min_age_column','max_age_column']
-    def __init__(self,*, min_age_column:Optional[str]='minAge', max_age_column:Optional[str]='maxAge'):
+    required_arguments = ['min_age_column','max_age_column','unit_name_column']
+    def __init__(self,*, unit_name_column:Optional[str]='name',min_age_column:Optional[str]='minAge', max_age_column:Optional[str]='maxAge'):
         """
         Initialiser for age based sorter
         """
         super().__init__()
+        self.unit_name_column = unit_name_column
         self.min_age_column = min_age_column
         self.max_age_column = max_age_column
         self.sorter_label = "SorterAgeBased"
@@ -204,9 +208,9 @@ def sort(self, units: pandas.DataFrame) -> list:
             sorted_units = sorted_units.sort_values(by=["meanAge"])
         logger.info("Stratigraphic order calculated using age based sorting")
         for _i, row in sorted_units.iterrows():
-            logger.info(f"{row['name']} - {row['minAge']} - {row['maxAge']}")
+            logger.info(f"{row[self.unit_name_column]} - {row[self.min_age_column]} - {row[self.max_age_column]}")
 
-        return list(sorted_units["name"])
+        return list(sorted_units[self.unit_name_column])
 
 
 class SorterAlpha(Sorter):
@@ -214,11 +218,12 @@ class SorterAlpha(Sorter):
     Sorter class which returns a sorted list of units based on the adjacency of units
     prioritising the units with lower number of contacting units
     """
-    required_arguments = ['contacts']
+    required_arguments = ['contacts', 'unit_name_column']
     def __init__(
         self,
         *,
         contacts: Optional[geopandas.GeoDataFrame] = None,
+        unit_name_column:Optional[str]='name',
     ):
         """
         Initialiser for adjacency based sorter
@@ -228,6 +233,7 @@ def __init__(
         """
         super().__init__()
         self.contacts = contacts
+        self.unit_name_column = unit_name_column
         self.sorter_label = "SorterAlpha"
         if 'UNITNAME_1' not in contacts.columns or 'UNITNAME_2' not in contacts.columns or 'length' not in contacts.columns:
             raise ValueError("contacts GeoDataFrame must contain 'UNITNAME_1', 'UNITNAME_2' and 'length' columns")
@@ -301,11 +307,14 @@ class SorterMaximiseContacts(Sorter):
     Sorter class which returns a sorted list of units based on the adjacency of units
     prioritising the maximum length of each contact
     """
-    required_arguments = ['contacts']
+    required_arguments = ['contacts', 'unit_name_column', 'unitname1_column', 'unitname2_column']
     def __init__(
         self,
         *,
         contacts: Optional[geopandas.GeoDataFrame] = None,
+        unit_name_column: str = 'name',
+        unitname1_column: str = 'UNITNAME_1',
+        unitname2_column: str = 'UNITNAME_2',
     ):
         """
         Initialiser for adjacency based sorter
@@ -320,8 +329,11 @@ def __init__(
         self.route = None
         self.directed_graph = None
         self.contacts = contacts
-        if 'UNITNAME_1' not in contacts.columns or 'UNITNAME_2' not in contacts.columns or 'length' not in contacts.columns:
-            raise ValueError("contacts GeoDataFrame must contain 'UNITNAME_1', 'UNITNAME_2' and 'length' columns")
+        self.unit_name_column = unit_name_column
+        self.unitname1_column = unitname1_column
+        self.unitname2_column = unitname2_column
+        if self.unitname1_column not in contacts.columns or self.unitname2_column not in contacts.columns or 'length' not in contacts.columns:
+            raise ValueError(f"contacts GeoDataFrame must contain '{self.unitname1_column}', '{self.unitname2_column}' and 'length' columns")
 
     def sort(self, units: pandas.DataFrame) -> list:
         """
@@ -339,13 +351,13 @@ def sort(self, units: pandas.DataFrame) -> list:
             raise ValueError("SorterMaximiseContacts requires 'contacts' argument")
         sorted_contacts = self.contacts.sort_values(by="length", ascending=False)
         self.graph = nx.Graph()
-        unit_names = list(units["name"].unique())
+        unit_names = list(units[self.unit_name_column].unique())
         for unit in unit_names:
             ## some units may not have any contacts e.g. if they are intrusives or sills. If we leave this then the
             ## sorter crashes
             if (
-                unit not in sorted_contacts['UNITNAME_1'].values
-                or unit not in sorted_contacts['UNITNAME_2'].values
+                unit not in sorted_contacts[self.unitname1_column].values
+                or unit not in sorted_contacts[self.unitname2_column].values
             ):
                 continue
             self.graph.add_node(unit, name=unit)
@@ -353,7 +365,7 @@ def sort(self, units: pandas.DataFrame) -> list:
         max_weight = max(list(sorted_contacts["length"])) + 1
         sorted_contacts['length'] /= max_weight
         for _, row in sorted_contacts.iterrows():
-            self.graph.add_edge(row["UNITNAME_1"], row["UNITNAME_2"], weight=(1 - row["length"]))
+            self.graph.add_edge(row[self.unitname1_column], row[self.unitname2_column], weight=(1 - row["length"]))
 
         self.route = nx_app.traveling_salesman_problem(self.graph)
         edge_list = list(nx.utils.pairwise(self.route))
@@ -384,10 +396,13 @@ class SorterObservationProjections(Sorter):
     Sorter class which returns a sorted list of units based on the adjacency of units
     using the direction of observations to predict which unit is adjacent to the current one
     """
-    required_arguments = ['contacts', 'geology_data', 'structure_data', 'dtm_data']
+    required_arguments = ['contacts', 'geology_data', 'structure_data', 'dtm_data', 'unit_name_column', 'unitname1_column', 'unitname2_column']
     def __init__(
         self,
         *,
+        unitname1_column: Optional[str] = 'UNITNAME_1',
+        unitname2_column: Optional[str] = 'UNITNAME_2',
+        unit_name_column: Optional[str] = 'name',
         contacts: Optional[geopandas.GeoDataFrame] = None,
         geology_data: Optional[geopandas.GeoDataFrame] = None,
         structure_data: Optional[geopandas.GeoDataFrame] = None,
@@ -409,9 +424,12 @@ def __init__(
         self.geology_data = geology_data
         self.structure_data = structure_data
         self.dtm_data = dtm_data
+        self.unit_name_column = unit_name_column
         self.sorter_label = "SorterObservationProjections"
         self.length = length
         self.lines = []
+        self.unit1name_column = unitname1_column
+        self.unit2name_column = unitname2_column
 
     def sort(self, units: pandas.DataFrame) -> list:
         """
@@ -543,13 +561,13 @@ def sort(self, units: pandas.DataFrame) -> list:
         g_undirected = g.to_undirected()
         for unit in unit_names:
             if len(list(g_undirected.neighbors(unit))) < 1:
-                mask1 = self.contacts["UNITNAME_1"] == unit
-                mask2 = self.contacts["UNITNAME_2"] == unit
+                mask1 = self.contacts[self.unit1name_column] == unit
+                mask2 = self.contacts[self.unit2name_column] == unit
                 for _, row in self.contacts[mask1 | mask2].iterrows():
-                    if unit == row["UNITNAME_1"]:
-                        g.add_edge(row["UNITNAME_2"], unit, weight=max_value * 10)
+                    if unit == row[self.unit1name_column]:
+                        g.add_edge(row[self.unit2name_column], unit, weight=max_value * 10)
                     else:
-                        g.add_edge(row["UNITNAME_1"], unit, weight=max_value * 10)
+                        g.add_edge(row[self.unit1name_column], unit, weight=max_value * 10)
 
         # Run travelling salesman using the observation evidence as weighting
         route = nx_app.traveling_salesman_problem(g.to_undirected())

From 9f3369fd91b41ccf8333bc869ee57a6b5ec19832 Mon Sep 17 00:00:00 2001
From: lachlangrose <lachlan.grose@monash.edu>
Date: Tue, 2 Dec 2025 10:49:09 +1100
Subject: [PATCH 11/31] fix: set age based sorter args

---
 map2loop/project.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/map2loop/project.py b/map2loop/project.py
index ec7f10f5..3cf355fd 100644
--- a/map2loop/project.py
+++ b/map2loop/project.py
@@ -614,6 +614,11 @@ def calculate_stratigraphic_order(self, take_best=False):
                 self.sorter.structure_data = self.map_data.get_map_data(Datatype.STRUCTURE)
             if hasattr(self.sorter, 'dtm_data') and self.sorter.dtm_data is None:
                 self.sorter.dtm_data = self.map_data.get_map_data(Datatype.DTM)
+            if hasattr(self.sorter, 'min_age_column') and self.sorter.min_age_column is None:
+                self.sorter.min_age_column = self.stratigraphic_column.get_min_age_column()
+            if hasattr(self.sorter, 'max_age_column') and self.sorter.max_age_column is None:
+                self.sorter.max_age_column = self.stratigraphic_column.get_max_age_column()
+            
             self.stratigraphic_column.column = self.sorter.sort(
                 self.stratigraphic_column.stratigraphicUnits,
             )

From 144a2249fb5ee34d94d8a47091b7c9f582c7ff28 Mon Sep 17 00:00:00 2001
From: lachlangrose <lachlan.grose@monash.edu>
Date: Tue, 2 Dec 2025 10:50:10 +1100
Subject: [PATCH 12/31] fix: project import for test

---
 tests/project/test_thickness_calculations.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/project/test_thickness_calculations.py b/tests/project/test_thickness_calculations.py
index 961e1a9c..7f601c84 100644
--- a/tests/project/test_thickness_calculations.py
+++ b/tests/project/test_thickness_calculations.py
@@ -5,7 +5,7 @@
 
 from map2loop.utils import load_map2loop_data
 from map2loop.thickness_calculator import InterpolatedStructure, StructuralPoint
-from map2loop import Project
+from map2loop.project import Project
 
 # 1. self.stratigraphic_column.stratigraphicUnits,
 st_units = pandas.DataFrame(

From ddbd0c4d7f5e457af53d019d6a0fef6609f204e4 Mon Sep 17 00:00:00 2001
From: lachlangrose <lachlan.grose@monash.edu>
Date: Tue, 2 Dec 2025 11:14:20 +1100
Subject: [PATCH 13/31] fix: updating requirements for alpha sorter plus add
 check for empty contacts

---
 map2loop/sorter.py | 28 ++++++++++++++++++----------
 1 file changed, 18 insertions(+), 10 deletions(-)

diff --git a/map2loop/sorter.py b/map2loop/sorter.py
index 738cfd25..b4aca705 100644
--- a/map2loop/sorter.py
+++ b/map2loop/sorter.py
@@ -9,6 +9,7 @@
 from osgeo import gdal
 from map2loop.utils import value_from_raster
 from .logging import getLogger
+import networkx as nx
 
 logger = getLogger(__name__)
 
@@ -218,12 +219,14 @@ class SorterAlpha(Sorter):
     Sorter class which returns a sorted list of units based on the adjacency of units
     prioritising the units with lower number of contacting units
     """
-    required_arguments = ['contacts', 'unit_name_column']
+    required_arguments = ['contacts', 'unit_name_column', 'unitname1_column', 'unitname2_column']
     def __init__(
         self,
         *,
         contacts: Optional[geopandas.GeoDataFrame] = None,
         unit_name_column:Optional[str]='name',
+        unitname1_column:Optional[str]='UNITNAME_1',
+        unitname2_column:Optional[str]='UNITNAME_2',
     ):
         """
         Initialiser for adjacency based sorter
@@ -235,9 +238,11 @@ def __init__(
         self.contacts = contacts
         self.unit_name_column = unit_name_column
         self.sorter_label = "SorterAlpha"
-        if 'UNITNAME_1' not in contacts.columns or 'UNITNAME_2' not in contacts.columns or 'length' not in contacts.columns:
-            raise ValueError("contacts GeoDataFrame must contain 'UNITNAME_1', 'UNITNAME_2' and 'length' columns")
-
+        self.unitname1_column = unitname1_column
+        self.unitname2_column = unitname2_column
+        if self.unitname1_column not in contacts.columns or self.unitname2_column not in contacts.columns or 'length' not in contacts.columns:
+            raise ValueError(f"contacts GeoDataFrame must contain '{self.unitname1_column}', '{self.unitname2_column}' and 'length' columns")
+        
     def sort(self, units: pandas.DataFrame) -> list:
         """
         Execute sorter method takes unit data and returns the sorted unit names based on this algorithm.
@@ -250,20 +255,22 @@ def sort(self, units: pandas.DataFrame) -> list:
         """
         if self.contacts is None:
             raise ValueError("contacts must be set (not None) before calling sort() in SorterAlpha.")
-        import networkx as nx
-        if self.contacts is None:
-            raise ValueError("SorterAlpha requires 'contacts' argument")
+        if len(self.contacts) == 0:
+            raise ValueError("contacts GeoDataFrame is empty in SorterAlpha.")
+        if 'length' not in self.contacts.columns:
+            self.contacts['length'] = self.contacts.geometry.length
+        self.contacts['length'] = self.contacts['length'].astype(float)
         sorted_contacts = self.contacts.sort_values(by="length", ascending=False)[
-            ["UNITNAME_1", "UNITNAME_2", "length"]
+            [self.unitname1_column, self.unitname2_column, "length"]
         ]
-        unit_names = list(units["name"].unique())
+        unit_names = list(units[self.unit_name_column].unique())
         graph = nx.Graph()
         for unit in unit_names:
             graph.add_node(unit, name=unit)
         max_weight = max(list(sorted_contacts["length"])) + 1
         for _, row in sorted_contacts.iterrows():
             graph.add_edge(
-                row["UNITNAME_1"], row["UNITNAME_2"], weight=int(max_weight - row["length"])
+                row[self.unitname1_column], row[self.unitname2_column], weight=int(max_weight - row["length"])
             )
 
         cnode = None
@@ -534,6 +541,7 @@ def sort(self, units: pandas.DataFrame) -> list:
         df = pandas.DataFrame(0, index=unit_names, columns=unit_names)
         for younger, older in ordered_unit_observations:
             df.loc[younger, older] += 1
+        print(df, df.max())
         max_value = max(df.max())
 
         # Using the older/younger matrix create a directed graph

From d71c0015851e33c051d1f7d3eb66009479c56f11 Mon Sep 17 00:00:00 2001
From: lachlangrose <lachlan.grose@monash.edu>
Date: Tue, 2 Dec 2025 11:15:38 +1100
Subject: [PATCH 14/31] fix: update maximise contacts to compute length and
 check for empty contacts

---
 map2loop/sorter.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/map2loop/sorter.py b/map2loop/sorter.py
index b4aca705..02343e3a 100644
--- a/map2loop/sorter.py
+++ b/map2loop/sorter.py
@@ -356,6 +356,11 @@ def sort(self, units: pandas.DataFrame) -> list:
         import networkx.algorithms.approximation as nx_app
         if self.contacts is None:
             raise ValueError("SorterMaximiseContacts requires 'contacts' argument")
+        if len(self.contacts) == 0:
+            raise ValueError("contacts GeoDataFrame is empty in SorterMaximiseContacts.")
+        if "length" not in self.contacts.columns:
+            self.contacts['length'] = self.contacts.geometry.length
+        self.contacts['length'] = self.contacts['length'].astype(float)
         sorted_contacts = self.contacts.sort_values(by="length", ascending=False)
         self.graph = nx.Graph()
         unit_names = list(units[self.unit_name_column].unique())

From 96972046130a0bb7733be1b871a8cea8d8af46ce Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Tue, 16 Dec 2025 10:52:42 +1100
Subject: [PATCH 15/31] fix: add option to set logger level and handler. Makes
 passing handler from qgis to python possible

---
 map2loop/logging.py | 50 +++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 50 insertions(+)

diff --git a/map2loop/logging.py b/map2loop/logging.py
index 816ac45e..b7e41b3f 100644
--- a/map2loop/logging.py
+++ b/map2loop/logging.py
@@ -66,3 +66,53 @@ def set_level(level: str):
         logger = logging.getLogger(name)
         logger.setLevel(level)
     logger.info(f"Logging level set to {level}")
+
+
+logger = getLogger(__name__)
+logger.info("Imported LoopStructural")
+
+
+def setLogging(level="info", handler=None):
+    """Set the logging parameters for log file or custom handler.
+
+    Parameters
+    ----------
+    level : str, optional
+        Logging level to set, by default "info"
+        Valid options: 'info', 'warning', 'error', 'debug'
+    handler : logging.Handler, optional
+        A logging handler to use instead of the default StreamHandler, by default None
+
+    Examples
+    --------
+    >>> from map2loop.logging import setLogging
+    >>> setLogging('debug')
+    >>> setLogging('info', logging.FileHandler('loop.log'))
+    """
+    levels = get_levels()
+    level_value = levels.get(level, logging.WARNING)
+
+    # Create default handler if none provided
+    if handler is None:
+        handler = logging.StreamHandler()
+
+    formatter = logging.Formatter(
+        "%(levelname)s: %(asctime)s: %(filename)s:%(lineno)d -- %(message)s"
+    )
+    handler.setFormatter(formatter)
+    handler.setLevel(level_value)
+
+    # Replace handlers in all known loggers
+    for name in loggers:
+        logger = logging.getLogger(name)
+        logger.handlers = []
+        logger.addHandler(handler)
+        logger.setLevel(level_value)
+
+    # Also apply to main module logger
+    main_logger = logging.getLogger(__name__)
+    main_logger.handlers = []
+    main_logger.addHandler(handler)
+    main_logger.setLevel(level_value)
+
+    main_logger.info(f"Set logging to {level}")

From c7435b517bbe98d7fb55ed4f11500308dfc455c0 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Tue, 16 Dec 2025 10:53:19 +1100
Subject: [PATCH 16/31] fix: updating incorrect type hints

---
 map2loop/utils/utility_functions.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/map2loop/utils/utility_functions.py b/map2loop/utils/utility_functions.py
index 082c615d..58c52721 100644
--- a/map2loop/utils/utility_functions.py
+++ b/map2loop/utils/utility_functions.py
@@ -184,7 +184,7 @@ def find_segment_strike_from_pt(
 
 @beartype.beartype
 def calculate_endpoints(
-    start_point: shapely.geometry.Point, azimuth_deg: float, distance: int, bbox: pandas.DataFrame
+    start_point: shapely.geometry.Point, azimuth_deg: float, distance: Union[float,int], bbox: pandas.DataFrame
 ) -> shapely.geometry.LineString:
     """
     Calculate the endpoints of a line segment given a start point, azimuth angle, distance, and bounding box.

From e8031c9a21a57d8eccd11706a71d212ff8904ca3 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Tue, 16 Dec 2025 10:53:54 +1100
Subject: [PATCH 17/31] fix: adding safeguarding and fixing merge

---
 map2loop/thickness_calculator.py | 30 ++++++++++++++++--------------
 1 file changed, 16 insertions(+), 14 deletions(-)

diff --git a/map2loop/thickness_calculator.py b/map2loop/thickness_calculator.py
index 83f2a290..9391462f 100644
--- a/map2loop/thickness_calculator.py
+++ b/map2loop/thickness_calculator.py
@@ -113,7 +113,7 @@ def _check_thickness_percentage_calculations(self, thicknesses: pandas.DataFrame
                 f"have a calculated thickness of -1. This may indicate that {self.thickness_calculator_label} "
                 f"is not suitable for this dataset."
             )
-            
+
 class ThicknessCalculatorAlpha(ThicknessCalculator):
     """
     ThicknessCalculator class which estimates unit thickness based on units, basal_contacts and stratigraphic order
@@ -452,8 +452,10 @@ def compute(
                 )
         
         # Combine all location_tracking DataFrames into a single DataFrame
-        combined_location_tracking = pandas.concat(_location_tracking, ignore_index=True)
-        
+        if _location_tracking and len(_location_tracking) > 0:
+            combined_location_tracking = pandas.concat(_location_tracking, ignore_index=True)
+        else:
+            combined_location_tracking = pandas.DataFrame()
         # Save the combined DataFrame as an attribute of the class
         self.location_tracking = combined_location_tracking
         
@@ -492,7 +494,6 @@ def __init__(
         self.strike_allowance = 30
         self.lines = None
 
-
     @beartype.beartype
     def compute(
         self,
@@ -552,7 +553,9 @@ def compute(
             crs=basal_contacts.crs,
         )
         # add unitname to the sampled structures
-        sampled_structures['unit_name'] = geopandas.sjoin(sampled_structures, geology)['UNITNAME']
+        sampled_structures['unit_name'] = geopandas.sjoin(
+            sampled_structures.drop(columns=['UNITNAME']), geology, how='inner', predicate='within'
+        )['UNITNAME']
 
         # remove nans from sampled structures
         # this happens when there are strati measurements within intrusions. If intrusions are removed from the geology map, unit_name will then return a nan
@@ -683,7 +686,7 @@ def compute(
             if not (b_s[0] < strike1 < b_s[1] and b_s[0] < strike2 < b_s[1]):
                 continue
 
-            #build the debug info
+            # build the debug info
             line = shapely.geometry.LineString([int_pt1, int_pt2])
             _lines.append(line)
             _dip.append(measurement['DIP'])  
@@ -694,17 +697,17 @@ def compute(
             # if length is higher than max_line_length, skip
             if self.max_line_length is not None and L > self.max_line_length:
                 continue
-            
+
             # calculate thickness
             thickness = L * math.sin(math.radians(measurement['DIP']))
 
             thicknesses.append(thickness)
             lis.append(litho_in)
-        
+
         # create the debug gdf
         self.lines = geopandas.GeoDataFrame(geometry=_lines, crs=basal_contacts.crs)
         self.lines["DIP"] = _dip
-        
+
         # create a DataFrame of the thicknesses median and standard deviation by lithology
         result = pandas.DataFrame({'unit': lis, 'thickness': thicknesses})
         result = result.groupby('unit')['thickness'].agg(['median', 'mean', 'std']).reset_index()
@@ -715,7 +718,7 @@ def compute(
         output_units['ThicknessMedian'] = numpy.full(len(output_units), numpy.nan)
         output_units['ThicknessMean'] = numpy.full(len(output_units), numpy.nan)
         output_units['ThicknessStdDev'] = numpy.full(len(output_units), numpy.nan)
-        
+
         # find which units have no thickness calculated
         names_not_in_result = units[~units['name'].isin(result['unit'])]['name'].to_list()
         # assign the thicknesses to the each unit
@@ -724,12 +727,11 @@ def compute(
             output_units.loc[idx, 'ThicknessMedian'] = unit['median']
             output_units.loc[idx, 'ThicknessMean'] = unit['mean']
             output_units.loc[idx, 'ThicknessStdDev'] = unit['std']
-       
+
         output_units["ThicknessMean"] = output_units["ThicknessMean"].fillna(-1)
         output_units["ThicknessMedian"] = output_units["ThicknessMedian"].fillna(-1)
         output_units["ThicknessStdDev"] = output_units["ThicknessStdDev"].fillna(-1)
-        
-        
+
         # handle the units that have no thickness
         for unit in names_not_in_result:
             # if no thickness has been calculated for the unit
@@ -760,7 +762,7 @@ def compute(
                 output_units.loc[output_units["name"] == unit, "ThicknessStdDev"] = -1
 
         self._check_thickness_percentage_calculations(output_units)
-        
+
         return output_units
 
 

From 7d581cfc2261b586657ce00c21b83941a28a2caa Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Tue, 16 Dec 2025 13:33:32 +1100
Subject: [PATCH 18/31] fix: save location tracking as a line geodataframe for
 qgis compatibility

---
 map2loop/thickness_calculator.py | 14 ++++++++++++--
 1 file changed, 12 insertions(+), 2 deletions(-)

diff --git a/map2loop/thickness_calculator.py b/map2loop/thickness_calculator.py
index 9391462f..2469e749 100644
--- a/map2loop/thickness_calculator.py
+++ b/map2loop/thickness_calculator.py
@@ -29,6 +29,7 @@
 import pandas
 import geopandas
 import shapely
+from shapely.geometry import LineString
 import math
 from osgeo import gdal
 from shapely.errors import UnsupportedGEOSVersionError
@@ -457,8 +458,17 @@ def compute(
         else:
             combined_location_tracking = pandas.DataFrame()
         # Save the combined DataFrame as an attribute of the class
-        self.location_tracking = combined_location_tracking
-        
+        # self.location_tracking = combined_location_tracking
+        combined_location_tracking['geometry'] = combined_location_tracking.apply(
+            lambda row: LineString([
+                (row['p1_x'], row['p1_y'], row['p1_z']),
+                (row['p2_x'], row['p2_y'], row['p2_z'])
+            ]),
+            axis=1
+        )
+
+        # Convert to GeoDataFrame and set CRS to match basal_contacts
+        self.location_tracking = geopandas.GeoDataFrame(combined_location_tracking, geometry='geometry', crs=basal_contacts.crs)
         # Create GeoDataFrame for lines
         self.lines = geopandas.GeoDataFrame(geometry=_lines, crs=basal_contacts.crs)
         self.lines['dip'] = _dips

From 33e2f66f5a4cbb45f438314d8b883b00ad52a013 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Tue, 16 Dec 2025 13:55:25 +1100
Subject: [PATCH 19/31] fix: only store geodataframe if layer is not empty

---
 map2loop/thickness_calculator.py | 48 +++++++++++++++++---------------
 1 file changed, 26 insertions(+), 22 deletions(-)

diff --git a/map2loop/thickness_calculator.py b/map2loop/thickness_calculator.py
index 2469e749..f75d00f5 100644
--- a/map2loop/thickness_calculator.py
+++ b/map2loop/thickness_calculator.py
@@ -239,7 +239,6 @@ def __init__(
         super().__init__(dtm_data, bounding_box, max_line_length, is_strike)
         self.thickness_calculator_label = "InterpolatedStructure"
         self.lines = None
-        
 
     @beartype.beartype
     def compute(
@@ -298,7 +297,7 @@ def compute(
         # set the crs of the contacts to the crs of the units
         contacts = contacts.set_crs(crs=basal_contacts.crs)
         if self.dtm_data is not None:
-        # get the elevation Z of the contacts
+            # get the elevation Z of the 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(
@@ -347,11 +346,11 @@ def compute(
         interpolated_orientations = interpolated_orientations[
             ["geometry", "dip", "UNITNAME"]
         ].copy()
-        
+
         _lines = []
         _dips = []
         _location_tracking = []
-        
+
         for i in range(0, len(stratigraphic_order) - 1):
             if (
                 stratigraphic_order[i] in basal_unit_list
@@ -373,21 +372,21 @@ def compute(
                     dip = interpolated_orientations.loc[
                         interpolated_orientations["UNITNAME"] == stratigraphic_order[i], "dip"
                     ].to_numpy()
-                    
+
                     _thickness = []
-                
+
                     for _, row in basal_contact.iterrows():
                         # find the shortest line between the basal contact points and top contact points
                         short_line = shapely.shortest_line(row.geometry, top_contact_geometry)
                         _lines.append(short_line[0])
-                        
-                        # check if the short line is 
+
+                        # check if the short line is
                         if self.max_line_length is not None and short_line.length > self.max_line_length:
                             continue
                         if self.dtm_data is not None:
                             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])
@@ -415,7 +414,7 @@ def compute(
                         _dips.append(_dip)
                         # calculate the true thickness t = L * sin(dip)
                         thickness = line_length * numpy.sin(_dip)
-                        
+
                         # add location tracking
                         location_tracking = pandas.DataFrame(
                             {
@@ -426,7 +425,7 @@ def compute(
                             }
                         )
                         _location_tracking.append(location_tracking)
-                        
+
                         # Average thickness along the shortest line
                         if all(numpy.isnan(thickness)):
                             pass
@@ -451,31 +450,36 @@ def compute(
                 logger.warning(
                     f"Thickness Calculator InterpolatedStructure: Cannot calculate thickness between {stratigraphic_order[i]} and {stratigraphic_order[i + 1]}\n"
                 )
-        
+
         # Combine all location_tracking DataFrames into a single DataFrame
         if _location_tracking and len(_location_tracking) > 0:
             combined_location_tracking = pandas.concat(_location_tracking, ignore_index=True)
+            combined_location_tracking['geometry'] = combined_location_tracking.apply(
+                lambda row: LineString(
+                    [
+                        (row['p1_x'], row['p1_y'], row['p1_z']),
+                        (row['p2_x'], row['p2_y'], row['p2_z']),
+                    ]
+                ),
+                axis=1,
+            )
         else:
             combined_location_tracking = pandas.DataFrame()
         # Save the combined DataFrame as an attribute of the class
         # self.location_tracking = combined_location_tracking
-        combined_location_tracking['geometry'] = combined_location_tracking.apply(
-            lambda row: LineString([
-                (row['p1_x'], row['p1_y'], row['p1_z']),
-                (row['p2_x'], row['p2_y'], row['p2_z'])
-            ]),
-            axis=1
-        )
 
         # Convert to GeoDataFrame and set CRS to match basal_contacts
-        self.location_tracking = geopandas.GeoDataFrame(combined_location_tracking, geometry='geometry', crs=basal_contacts.crs)
+        if not combined_location_tracking.empty:
+            self.location_tracking = geopandas.GeoDataFrame(combined_location_tracking, geometry='geometry', crs=basal_contacts.crs)
+        else:
+            self.location_tracking = geopandas.GeoDataFrame(columns=['p1_x', 'p1_y', 'p1_z', 'p2_x', 'p2_y', 'p2_z', 'thickness', 'unit', 'geometry'], crs=basal_contacts.crs)
         # Create GeoDataFrame for lines
         self.lines = geopandas.GeoDataFrame(geometry=_lines, crs=basal_contacts.crs)
         self.lines['dip'] = _dips
-        
+
         # Check thickness calculation
         self._check_thickness_percentage_calculations(thicknesses)
-        
+
         return thicknesses
 
 class StructuralPoint(ThicknessCalculator):

From 4f4d0a5f2bfe96dff2bf48dd766d53b430280c57 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Wed, 17 Dec 2025 10:13:16 +1100
Subject: [PATCH 20/31] fix: updating logic

---
 map2loop/thickness_calculator.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/map2loop/thickness_calculator.py b/map2loop/thickness_calculator.py
index f75d00f5..a7ef7cf4 100644
--- a/map2loop/thickness_calculator.py
+++ b/map2loop/thickness_calculator.py
@@ -566,9 +566,11 @@ def compute(
             geometry=geopandas.points_from_xy(sampled_structures.X, sampled_structures.Y),
             crs=basal_contacts.crs,
         )
+        if 'UNITNAME' in sampled_structures.columns:
+            sampled_structures = sampled_structures.drop(columns=['UNITNAME'])
         # add unitname to the sampled structures
         sampled_structures['unit_name'] = geopandas.sjoin(
-            sampled_structures.drop(columns=['UNITNAME']), geology, how='inner', predicate='within'
+            sampled_structures, geology, how='inner', predicate='within'
         )['UNITNAME']
 
         # remove nans from sampled structures

From b8c71f02443a91aa8ad1f028f4f42c33de365972 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Thu, 18 Dec 2025 08:26:31 +1100
Subject: [PATCH 21/31] fix: beartype issues

---
 map2loop/thickness_calculator.py    | 2 +-
 map2loop/utils/utility_functions.py | 6 +++---
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/map2loop/thickness_calculator.py b/map2loop/thickness_calculator.py
index a7ef7cf4..82d7b607 100644
--- a/map2loop/thickness_calculator.py
+++ b/map2loop/thickness_calculator.py
@@ -630,7 +630,7 @@ def compute(
             # draw orthogonal line to the strike (default value 10Km), and clip it by the bounding box of the lithology
             if self.max_line_length is None:
                 self.max_line_length = 10000
-            B = calculate_endpoints(measurement_pt, strike, self.max_line_length, bbox_poly)
+            B = calculate_endpoints(measurement_pt, float(strike), float(self.max_line_length), bbox_poly)
             b = geopandas.GeoDataFrame({'geometry': [B]}).set_crs(basal_contacts.crs)
 
             # find all intersections
diff --git a/map2loop/utils/utility_functions.py b/map2loop/utils/utility_functions.py
index 58c52721..ea737a3c 100644
--- a/map2loop/utils/utility_functions.py
+++ b/map2loop/utils/utility_functions.py
@@ -184,8 +184,8 @@ def find_segment_strike_from_pt(
 
 @beartype.beartype
 def calculate_endpoints(
-    start_point: shapely.geometry.Point, azimuth_deg: float, distance: Union[float,int], bbox: pandas.DataFrame
-) -> shapely.geometry.LineString:
+    start_point: shapely.geometry.Point, azimuth_deg: Union[float,int], distance: Union[float,int], bbox: pandas.DataFrame
+) -> Union[shapely.geometry.LineString,shapely.geometry.GeometryCollection]:
     """
     Calculate the endpoints of a line segment given a start point, azimuth angle, distance, and bounding box.
 
@@ -220,7 +220,7 @@ def calculate_endpoints(
     line = shapely.geometry.LineString([left_endpoint, right_endpoint])
 
     new_line = shapely.ops.clip_by_rect(line, minx, miny, maxx, maxy)
-
+    print(type(new_line))
     return new_line
 
 

From 1e4882d44644cc72fc44bae3a54e5d98f41db20c Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Thu, 18 Dec 2025 08:36:05 +1100
Subject: [PATCH 22/31] fix: add ability to specify fault field to topoloy

---
 map2loop/topology.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/map2loop/topology.py b/map2loop/topology.py
index d7a90dc4..88b4e643 100644
--- a/map2loop/topology.py
+++ b/map2loop/topology.py
@@ -36,7 +36,8 @@ class Topology:
     def __init__(
         self, geology_data: geopandas.GeoDataFrame, 
         fault_data: Optional[geopandas.GeoDataFrame] = None,  
-        verbose_level: VerboseLevel = VerboseLevel.NONE
+        verbose_level: VerboseLevel = VerboseLevel.NONE,
+        id_field: str = 'ID',
     ):
         """
         The initialiser for the map2model wrapper
@@ -55,7 +56,7 @@ def __init__(
         self.fault_data = fault_data
         self.verbose_level = verbose_level
         self.buffer_radius = 500
-
+        self.fault_id_field = id_field
     @property
     def fault_fault_relationships(self):
         if self._fault_fault_relationships is None:

From c44913b2dfae93d50865cd2247191faf067b6a95 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Thu, 18 Dec 2025 17:07:09 +1100
Subject: [PATCH 23/31] fix: add call methods to allow for generic debug
 exporting

---
 map2loop/sampler.py              | 3 ++-
 map2loop/sorter.py               | 4 +++-
 map2loop/thickness_calculator.py | 2 ++
 3 files changed, 7 insertions(+), 2 deletions(-)

diff --git a/map2loop/sampler.py b/map2loop/sampler.py
index b4c7835c..0544f2f5 100644
--- a/map2loop/sampler.py
+++ b/map2loop/sampler.py
@@ -52,7 +52,8 @@ def sample(
             pandas.DataFrame: data frame containing samples
         """
         pass
-
+    def __call__(self, *args):
+        return self.sample(*args)
 
 class SamplerDecimator(Sampler):
     """
diff --git a/map2loop/sorter.py b/map2loop/sorter.py
index 02343e3a..730b2f41 100644
--- a/map2loop/sorter.py
+++ b/map2loop/sorter.py
@@ -59,7 +59,9 @@ def sort(self, units: pandas.DataFrame) -> list:
             list: sorted list of unit names
         """
         pass
-
+    
+    def __call__(self, *args):
+        return self.sort(*args)
 
 class SorterUseNetworkX(Sorter):
     """
diff --git a/map2loop/thickness_calculator.py b/map2loop/thickness_calculator.py
index 82d7b607..1377a520 100644
--- a/map2loop/thickness_calculator.py
+++ b/map2loop/thickness_calculator.py
@@ -114,6 +114,8 @@ def _check_thickness_percentage_calculations(self, thicknesses: pandas.DataFrame
                 f"have a calculated thickness of -1. This may indicate that {self.thickness_calculator_label} "
                 f"is not suitable for this dataset."
             )
+    def __call__(self, *args):
+        return self.compute(*args)
 
 class ThicknessCalculatorAlpha(ThicknessCalculator):
     """

From a5ed09594d8745d3399de42c6631fc4e383562cd Mon Sep 17 00:00:00 2001
From: lachlangrose <lachlan.grose@monash.edu>
Date: Mon, 19 Jan 2026 10:57:24 +1100
Subject: [PATCH 24/31] fix: replace args for kwargs

---
 map2loop/sampler.py              | 4 ++--
 map2loop/sorter.py               | 4 ++--
 map2loop/thickness_calculator.py | 4 ++--
 3 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/map2loop/sampler.py b/map2loop/sampler.py
index 0544f2f5..6aad3550 100644
--- a/map2loop/sampler.py
+++ b/map2loop/sampler.py
@@ -52,8 +52,8 @@ def sample(
             pandas.DataFrame: data frame containing samples
         """
         pass
-    def __call__(self, *args):
-        return self.sample(*args)
+    def __call__(self, **kwargs):
+        return self.sample(**kwargs)
 
 class SamplerDecimator(Sampler):
     """
diff --git a/map2loop/sorter.py b/map2loop/sorter.py
index 730b2f41..4e162119 100644
--- a/map2loop/sorter.py
+++ b/map2loop/sorter.py
@@ -60,8 +60,8 @@ def sort(self, units: pandas.DataFrame) -> list:
         """
         pass
     
-    def __call__(self, *args):
-        return self.sort(*args)
+    def __call__(self, **kwargs):
+        return self.sort(**kwargs)
 
 class SorterUseNetworkX(Sorter):
     """
diff --git a/map2loop/thickness_calculator.py b/map2loop/thickness_calculator.py
index 1377a520..f5585dd4 100644
--- a/map2loop/thickness_calculator.py
+++ b/map2loop/thickness_calculator.py
@@ -114,8 +114,8 @@ def _check_thickness_percentage_calculations(self, thicknesses: pandas.DataFrame
                 f"have a calculated thickness of -1. This may indicate that {self.thickness_calculator_label} "
                 f"is not suitable for this dataset."
             )
-    def __call__(self, *args):
-        return self.compute(*args)
+    def __call__(self, **kwargs):
+        return self.compute(**kwargs)
 
 class ThicknessCalculatorAlpha(ThicknessCalculator):
     """

From 651a9e619385e5754556e231818dbc9140aea545 Mon Sep 17 00:00:00 2001
From: lachlangrose <lachlan.grose@monash.edu>
Date: Mon, 19 Jan 2026 11:08:42 +1100
Subject: [PATCH 25/31] fix: applying copilot suggestions

---
 map2loop/logging.py                 | 2 +-
 map2loop/utils/utility_functions.py | 1 -
 2 files changed, 1 insertion(+), 2 deletions(-)

diff --git a/map2loop/logging.py b/map2loop/logging.py
index b7e41b3f..c0fe217c 100644
--- a/map2loop/logging.py
+++ b/map2loop/logging.py
@@ -69,7 +69,7 @@ def set_level(level: str):
 
 
 logger = getLogger(__name__)
-logger.info("Imported LoopStructural")
+logger.info("Imported map2loop.logging module")
 
 
 def setLogging(level="info", handler=None):
diff --git a/map2loop/utils/utility_functions.py b/map2loop/utils/utility_functions.py
index ea737a3c..edc87468 100644
--- a/map2loop/utils/utility_functions.py
+++ b/map2loop/utils/utility_functions.py
@@ -220,7 +220,6 @@ def calculate_endpoints(
     line = shapely.geometry.LineString([left_endpoint, right_endpoint])
 
     new_line = shapely.ops.clip_by_rect(line, minx, miny, maxx, maxy)
-    print(type(new_line))
     return new_line
 
 

From ac280725f20409e8e2cab2cb19b2d3fe7ac12286 Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@monash.edu>
Date: Mon, 19 Jan 2026 11:01:23 +1100
Subject: [PATCH 26/31] Enable auto-update for conda in workflow

---
 .github/workflows/linting_and_testing.yml | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/.github/workflows/linting_and_testing.yml b/.github/workflows/linting_and_testing.yml
index 0259ccb5..1af945f6 100644
--- a/.github/workflows/linting_and_testing.yml
+++ b/.github/workflows/linting_and_testing.yml
@@ -46,9 +46,8 @@ jobs:
       - uses: actions/checkout@v4
       - uses: conda-incubator/setup-miniconda@v3
         with:
+            auto-update-conda: true
             python-version: ${{ matrix.python-version }}
-            conda-remove-defaults: "true"
-      
       - name: Install dependencies for windows python 3.9 and 3.10
         if: ${{ matrix.os == 'windows-latest' && (matrix.python-version == '3.9' || matrix.python-version == '3.10') }}
         run: |
@@ -77,4 +76,4 @@ jobs:
 
       - name: Run tests
         run: |
-            conda run -n test  pytest -v
\ No newline at end of file
+            conda run -n test  pytest -v

From 7cc4b099fe702d255f858770f0e3e4c247ea580a Mon Sep 17 00:00:00 2001
From: lachlangrose <lachlan.grose@monash.edu>
Date: Mon, 19 Jan 2026 11:11:09 +1100
Subject: [PATCH 27/31] style: formatting

---
 map2loop/sorter.py                  | 135 +++++++++++++++++++---------
 map2loop/thickness_calculator.py    |   1 -
 map2loop/utils/utility_functions.py |  45 ++++++----
 3 files changed, 117 insertions(+), 64 deletions(-)

diff --git a/map2loop/sorter.py b/map2loop/sorter.py
index 4e162119..123d357c 100644
--- a/map2loop/sorter.py
+++ b/map2loop/sorter.py
@@ -22,11 +22,10 @@ class Sorter(ABC):
         ABC (ABC): Derived from Abstract Base Class
     """
 
-    def __init__(
-        self):
+    def __init__(self):
         """
         Initialiser for Sorter
-        
+
         Args:
             unit_relationships (pandas.DataFrame): the relationships between units (columns must contain ["Index1", "Unitname1", "Index2", "Unitname2"])
             contacts (pandas.DataFrame): unit contacts with length of the contacts in metres
@@ -35,7 +34,6 @@ def __init__(
             dtm_data (gdal.Dataset): the dtm data
         """
         self.sorter_label = "SorterBaseClass"
-        
 
     def type(self):
         """
@@ -59,29 +57,28 @@ def sort(self, units: pandas.DataFrame) -> list:
             list: sorted list of unit names
         """
         pass
-    
+
     def __call__(self, **kwargs):
         return self.sort(**kwargs)
 
+
 class SorterUseNetworkX(Sorter):
     """
     Sorter class which returns a sorted list of units based on the unit relationships using a topological graph sorting algorithm
     """
-    required_arguments: List[str] = [
-        'geology_data',
-        'unit_name_column'
-    ]
+
+    required_arguments: List[str] = ['geology_data', 'unit_name_column']
 
     def __init__(
         self,
         *,
-        unit_name_column:Optional[str]='name',
+        unit_name_column: Optional[str] = 'name',
         unit_relationships: Optional[pandas.DataFrame] = None,
         geology_data: Optional[geopandas.GeoDataFrame] = None,
     ):
         """
         Initialiser for networkx graph sorter
-        
+
         Args:
             unit_relationships (pandas.DataFrame): the relationships between units
         """
@@ -94,6 +91,7 @@ def __init__(
             self.unit_relationships = unit_relationships
         else:
             self.unit_relationships = None
+
     def set_geology_data(self, geology_data: geopandas.GeoDataFrame):
         """
         Set geology data and calculate topology and unit relationships
@@ -102,19 +100,20 @@ def set_geology_data(self, geology_data: geopandas.GeoDataFrame):
             geology_data (geopandas.GeoDataFrame): the geology data
         """
         self._calculate_topology(geology_data)
+
     def _calculate_topology(self, geology_data: geopandas.GeoDataFrame):
         if geology_data is None:
             raise ValueError("geology_data is required")
 
         if isinstance(geology_data, geopandas.GeoDataFrame) is False:
             raise TypeError("geology_data must be a geopandas.GeoDataFrame")
-        
+
         if 'UNITNAME' not in geology_data.columns:
             raise ValueError("geology_data must contain 'UNITNAME' column")
-        
+
         self.topology = Topology(geology_data=geology_data)
         self.unit_relationships = self.topology.get_unit_unit_relationships()
-    
+
     @beartype.beartype
     def sort(self, units: pandas.DataFrame) -> list:
         """
@@ -127,6 +126,7 @@ def sort(self, units: pandas.DataFrame) -> list:
             list: the sorted unit names
         """
         import networkx as nx
+
         if self.unit_relationships is None:
             raise ValueError("SorterUseNetworkX requires 'unit_relationships' argument")
         graph = nx.DiGraph()
@@ -156,24 +156,26 @@ def sort(self, units: pandas.DataFrame) -> list:
 
 class SorterUseHint(SorterUseNetworkX):
     required_arguments: List[str] = ['unit_relationships']
-    def __init__(
-        self,
-        *,
-        geology_data: Optional[geopandas.GeoDataFrame] = None,
-    ):
-        logger.warning(
-            "SorterUseHint is deprecated in v3.2. Using SorterUseNetworkX instead"
-        )
+
+    def __init__(self, *, geology_data: Optional[geopandas.GeoDataFrame] = None):
+        logger.warning("SorterUseHint is deprecated in v3.2. Using SorterUseNetworkX instead")
         super().__init__(geology_data=geology_data)
-    
 
 
 class SorterAgeBased(Sorter):
     """
     Sorter class which returns a sorted list of units based on the min and max ages of the units
     """
-    required_arguments = ['min_age_column','max_age_column','unit_name_column']
-    def __init__(self,*, unit_name_column:Optional[str]='name',min_age_column:Optional[str]='minAge', max_age_column:Optional[str]='maxAge'):
+
+    required_arguments = ['min_age_column', 'max_age_column', 'unit_name_column']
+
+    def __init__(
+        self,
+        *,
+        unit_name_column: Optional[str] = 'name',
+        min_age_column: Optional[str] = 'minAge',
+        max_age_column: Optional[str] = 'maxAge',
+    ):
         """
         Initialiser for age based sorter
         """
@@ -202,16 +204,22 @@ def sort(self, units: pandas.DataFrame) -> list:
                 lambda row: (row[self.min_age_column] + row[self.max_age_column]) / 2.0, axis=1
             )
         else:
-            logger.error(f"Columns {self.min_age_column} and {self.max_age_column} must be present in units DataFrame")
+            logger.error(
+                f"Columns {self.min_age_column} and {self.max_age_column} must be present in units DataFrame"
+            )
             logger.error(f"Available columns are: {units.columns.tolist()}")
-            raise ValueError(f"Columns {self.min_age_column} and {self.max_age_column} must be present in units DataFrame")
+            raise ValueError(
+                f"Columns {self.min_age_column} and {self.max_age_column} must be present in units DataFrame"
+            )
         if "group" in units.columns:
             sorted_units = sorted_units.sort_values(by=["group", "meanAge"])
         else:
             sorted_units = sorted_units.sort_values(by=["meanAge"])
         logger.info("Stratigraphic order calculated using age based sorting")
         for _i, row in sorted_units.iterrows():
-            logger.info(f"{row[self.unit_name_column]} - {row[self.min_age_column]} - {row[self.max_age_column]}")
+            logger.info(
+                f"{row[self.unit_name_column]} - {row[self.min_age_column]} - {row[self.max_age_column]}"
+            )
 
         return list(sorted_units[self.unit_name_column])
 
@@ -221,18 +229,20 @@ class SorterAlpha(Sorter):
     Sorter class which returns a sorted list of units based on the adjacency of units
     prioritising the units with lower number of contacting units
     """
+
     required_arguments = ['contacts', 'unit_name_column', 'unitname1_column', 'unitname2_column']
+
     def __init__(
         self,
         *,
         contacts: Optional[geopandas.GeoDataFrame] = None,
-        unit_name_column:Optional[str]='name',
-        unitname1_column:Optional[str]='UNITNAME_1',
-        unitname2_column:Optional[str]='UNITNAME_2',
+        unit_name_column: Optional[str] = 'name',
+        unitname1_column: Optional[str] = 'UNITNAME_1',
+        unitname2_column: Optional[str] = 'UNITNAME_2',
     ):
         """
         Initialiser for adjacency based sorter
-        
+
         Args:
             contacts (geopandas.GeoDataFrame): unit contacts with length of the contacts in metres
         """
@@ -242,9 +252,15 @@ def __init__(
         self.sorter_label = "SorterAlpha"
         self.unitname1_column = unitname1_column
         self.unitname2_column = unitname2_column
-        if self.unitname1_column not in contacts.columns or self.unitname2_column not in contacts.columns or 'length' not in contacts.columns:
-            raise ValueError(f"contacts GeoDataFrame must contain '{self.unitname1_column}', '{self.unitname2_column}' and 'length' columns")
-        
+        if (
+            self.unitname1_column not in contacts.columns
+            or self.unitname2_column not in contacts.columns
+            or 'length' not in contacts.columns
+        ):
+            raise ValueError(
+                f"contacts GeoDataFrame must contain '{self.unitname1_column}', '{self.unitname2_column}' and 'length' columns"
+            )
+
     def sort(self, units: pandas.DataFrame) -> list:
         """
         Execute sorter method takes unit data and returns the sorted unit names based on this algorithm.
@@ -256,7 +272,9 @@ def sort(self, units: pandas.DataFrame) -> list:
             list: the sorted unit names
         """
         if self.contacts is None:
-            raise ValueError("contacts must be set (not None) before calling sort() in SorterAlpha.")
+            raise ValueError(
+                "contacts must be set (not None) before calling sort() in SorterAlpha."
+            )
         if len(self.contacts) == 0:
             raise ValueError("contacts GeoDataFrame is empty in SorterAlpha.")
         if 'length' not in self.contacts.columns:
@@ -272,7 +290,9 @@ def sort(self, units: pandas.DataFrame) -> list:
         max_weight = max(list(sorted_contacts["length"])) + 1
         for _, row in sorted_contacts.iterrows():
             graph.add_edge(
-                row[self.unitname1_column], row[self.unitname2_column], weight=int(max_weight - row["length"])
+                row[self.unitname1_column],
+                row[self.unitname2_column],
+                weight=int(max_weight - row["length"]),
             )
 
         cnode = None
@@ -316,7 +336,9 @@ class SorterMaximiseContacts(Sorter):
     Sorter class which returns a sorted list of units based on the adjacency of units
     prioritising the maximum length of each contact
     """
+
     required_arguments = ['contacts', 'unit_name_column', 'unitname1_column', 'unitname2_column']
+
     def __init__(
         self,
         *,
@@ -327,7 +349,7 @@ def __init__(
     ):
         """
         Initialiser for adjacency based sorter
-        
+
         Args:
             contacts (pandas.DataFrame): unit contacts with length of the contacts in metres
         """
@@ -341,8 +363,14 @@ def __init__(
         self.unit_name_column = unit_name_column
         self.unitname1_column = unitname1_column
         self.unitname2_column = unitname2_column
-        if self.unitname1_column not in contacts.columns or self.unitname2_column not in contacts.columns or 'length' not in contacts.columns:
-            raise ValueError(f"contacts GeoDataFrame must contain '{self.unitname1_column}', '{self.unitname2_column}' and 'length' columns")
+        if (
+            self.unitname1_column not in contacts.columns
+            or self.unitname2_column not in contacts.columns
+            or 'length' not in contacts.columns
+        ):
+            raise ValueError(
+                f"contacts GeoDataFrame must contain '{self.unitname1_column}', '{self.unitname2_column}' and 'length' columns"
+            )
 
     def sort(self, units: pandas.DataFrame) -> list:
         """
@@ -356,6 +384,7 @@ def sort(self, units: pandas.DataFrame) -> list:
         """
         import networkx as nx
         import networkx.algorithms.approximation as nx_app
+
         if self.contacts is None:
             raise ValueError("SorterMaximiseContacts requires 'contacts' argument")
         if len(self.contacts) == 0:
@@ -379,7 +408,9 @@ def sort(self, units: pandas.DataFrame) -> list:
         max_weight = max(list(sorted_contacts["length"])) + 1
         sorted_contacts['length'] /= max_weight
         for _, row in sorted_contacts.iterrows():
-            self.graph.add_edge(row[self.unitname1_column], row[self.unitname2_column], weight=(1 - row["length"]))
+            self.graph.add_edge(
+                row[self.unitname1_column], row[self.unitname2_column], weight=(1 - row["length"])
+            )
 
         self.route = nx_app.traveling_salesman_problem(self.graph)
         edge_list = list(nx.utils.pairwise(self.route))
@@ -410,7 +441,17 @@ class SorterObservationProjections(Sorter):
     Sorter class which returns a sorted list of units based on the adjacency of units
     using the direction of observations to predict which unit is adjacent to the current one
     """
-    required_arguments = ['contacts', 'geology_data', 'structure_data', 'dtm_data', 'unit_name_column', 'unitname1_column', 'unitname2_column']
+
+    required_arguments = [
+        'contacts',
+        'geology_data',
+        'structure_data',
+        'dtm_data',
+        'unit_name_column',
+        'unitname1_column',
+        'unitname2_column',
+    ]
+
     def __init__(
         self,
         *,
@@ -421,7 +462,7 @@ def __init__(
         geology_data: Optional[geopandas.GeoDataFrame] = None,
         structure_data: Optional[geopandas.GeoDataFrame] = None,
         dtm_data: Optional[gdal.Dataset] = None,
-        length: Union[float, int] = 1000
+        length: Union[float, int] = 1000,
     ):
         """
         Initialiser for adjacency based sorter
@@ -458,6 +499,7 @@ def sort(self, units: pandas.DataFrame) -> list:
         import networkx as nx
         import networkx.algorithms.approximation as nx_app
         from shapely.geometry import LineString, Point
+
         if self.contacts is None:
             raise ValueError("SorterObservationProjections requires 'contacts' argument")
         if self.geology_data is None:
@@ -525,7 +567,12 @@ def sort(self, units: pandas.DataFrame) -> list:
                     # Get heights for intersection point and start of ray
                     height = value_from_raster(inv_geotransform, dtm_array, start.x, start.y)
                     first_intersect_point = Point(start.x, start.y, height)
-                    height = value_from_raster(inv_geotransform, dtm_array, second_intersect_point.x, second_intersect_point.y)
+                    height = value_from_raster(
+                        inv_geotransform,
+                        dtm_array,
+                        second_intersect_point.x,
+                        second_intersect_point.y,
+                    )
                     second_intersect_point = Point(second_intersect_point.x, start.y, height)
 
                     # Check vertical difference between points and compare to projected dip angle
diff --git a/map2loop/thickness_calculator.py b/map2loop/thickness_calculator.py
index f5585dd4..5f6e3b57 100644
--- a/map2loop/thickness_calculator.py
+++ b/map2loop/thickness_calculator.py
@@ -468,7 +468,6 @@ def compute(
         else:
             combined_location_tracking = pandas.DataFrame()
         # Save the combined DataFrame as an attribute of the class
-        # self.location_tracking = combined_location_tracking
 
         # Convert to GeoDataFrame and set CRS to match basal_contacts
         if not combined_location_tracking.empty:
diff --git a/map2loop/utils/utility_functions.py b/map2loop/utils/utility_functions.py
index edc87468..42287d28 100644
--- a/map2loop/utils/utility_functions.py
+++ b/map2loop/utils/utility_functions.py
@@ -10,6 +10,7 @@
 from osgeo import gdal
 
 from ..logging import getLogger
+
 logger = getLogger(__name__)
 
 
@@ -184,8 +185,11 @@ def find_segment_strike_from_pt(
 
 @beartype.beartype
 def calculate_endpoints(
-    start_point: shapely.geometry.Point, azimuth_deg: Union[float,int], distance: Union[float,int], bbox: pandas.DataFrame
-) -> Union[shapely.geometry.LineString,shapely.geometry.GeometryCollection]:
+    start_point: shapely.geometry.Point,
+    azimuth_deg: Union[float, int],
+    distance: Union[float, int],
+    bbox: pandas.DataFrame,
+) -> Union[shapely.geometry.LineString, shapely.geometry.GeometryCollection]:
     """
     Calculate the endpoints of a line segment given a start point, azimuth angle, distance, and bounding box.
 
@@ -225,7 +229,7 @@ def calculate_endpoints(
 
 @beartype.beartype
 def multiline_to_line(
-    geometry: Union[shapely.geometry.LineString, shapely.geometry.MultiLineString]
+    geometry: Union[shapely.geometry.LineString, shapely.geometry.MultiLineString],
 ) -> shapely.geometry.LineString:
     """
     Converts a multiline geometry to a single line geometry.
@@ -382,7 +386,6 @@ def hex_to_rgb(hex_color: str) -> tuple:
 def calculate_minimum_fault_length(
     bbox: dict[str, Union[int, float]], area_percentage: float
 ) -> float:
-
     """
     Calculate the minimum fault length based on the map bounding box and a given area percentage.
 
@@ -440,11 +443,10 @@ def read_hjson_with_json(file_path: str) -> dict:
     except json.JSONDecodeError as e:
         raise ValueError(f"Failed to decode preprocessed HJSON as JSON: {e}") from e
 
+
 @beartype.beartype
 def update_from_legacy_file(
-    filename: str,
-    json_save_path: Optional[str] = None,
-    lower: bool = False
+    filename: str, json_save_path: Optional[str] = None, lower: bool = False
 ) -> Optional[Dict[str, Dict]]:
     """
     Update the config dictionary from the provided old version dictionary
@@ -452,18 +454,19 @@ def update_from_legacy_file(
         filename (str): the path to the legacy file
         json_save_path (str, optional): the path to save the updated json file. Defaults to None.
         lower (bool, optional): whether to convert all strings to lowercase. Defaults to False.
-        
+
     Returns:
         Dict[Dict]: the updated config dictionary
-    
+
     Example:
         from map2loop.utils import update_from_legacy_file
         update_from_legacy_file(filename=r"./source_data/example.hjson")
     """
     # only import config if needed
     from .config import Config
+
     file_map = Config()
-    
+
     code_mapping = {
         "otype": (file_map.structure_config, "orientation_type"),
         "dd": (file_map.structure_config, "dipdir_column"),
@@ -505,7 +508,7 @@ def update_from_legacy_file(
     except Exception as e:
         logger.error(f"Error reading file {filename}: {e}")
         return
-    #map the keys 
+    # map the keys
     file_map = file_map.to_dict()
     for legacy_key, new_mapping in code_mapping.items():
         if legacy_key in parsed_data:
@@ -514,21 +517,22 @@ def update_from_legacy_file(
             if lower and isinstance(value, str):
                 value = value.lower()
             section[new_key] = value
-    
+
     if "o" in parsed_data:
         object_id_value = parsed_data["o"]
         if lower and isinstance(object_id_value, str):
             object_id_value = object_id_value.lower()
         file_map['structure']["objectid_column"] = object_id_value
         file_map['geology']["objectid_column"] = object_id_value
-        file_map['fold']["objectid_column"] = object_id_value 
-        
+        file_map['fold']["objectid_column"] = object_id_value
+
     if json_save_path is not None:
         with open(json_save_path, "w") as f:
             json.dump(parsed_data, f, indent=4)
-    
+
     return file_map
 
+
 @beartype.beartype
 def value_from_raster(inv_geotransform, data, x: float, y: float):
     """
@@ -556,6 +560,7 @@ def value_from_raster(inv_geotransform, data, x: float, y: float):
     py = min(py, data.shape[1] - 1)
     return data[px][py]
 
+
 @beartype.beartype
 def set_z_values_from_raster_df(dtm_data: gdal.Dataset, df: pandas.DataFrame):
     """
@@ -573,7 +578,7 @@ def set_z_values_from_raster_df(dtm_data: gdal.Dataset, df: pandas.DataFrame):
     if len(df) <= 0:
         df["Z"] = []
         return df
-    
+
     if dtm_data is None:
         logger.warning("Cannot get value from data as data is not loaded")
         return None
@@ -582,11 +587,11 @@ def set_z_values_from_raster_df(dtm_data: gdal.Dataset, df: pandas.DataFrame):
     data_array = numpy.array(dtm_data.GetRasterBand(1).ReadAsArray().T)
 
     df["Z"] = df.apply(
-        lambda row: value_from_raster(inv_geotransform, data_array, row["X"], row["Y"]),
-        axis=1,
+        lambda row: value_from_raster(inv_geotransform, data_array, row["X"], row["Y"]), axis=1
     )
 
     return df
+<<<<<<< HEAD
 
 def clean_line_geometry(geometry: shapely.geometry.base.BaseGeometry):
     """Clean and normalize a Shapely geometry to a single LineString or None.
@@ -795,4 +800,6 @@ def segment_measure_range(parent_line: shapely.geometry.LineString, segment: sha
     end_measure = parent_line.project(end_point)
     if end_measure < start_measure:
         start_measure, end_measure = end_measure, start_measure
-    return start_measure, end_measure
\ No newline at end of file
+    return start_measure, end_measure
+=======
+>>>>>>> 1febaa9 (style: formatting)

From 4aff960243e8871c69a64604ee7b0f15fd4add25 Mon Sep 17 00:00:00 2001
From: "github-actions[bot]"
 <41898282+github-actions[bot]@users.noreply.github.com>
Date: Tue, 20 Jan 2026 23:41:47 +0000
Subject: [PATCH 28/31] chore(master): release 3.2.10

---
 .release-please-manifest.json |  2 +-
 CHANGELOG.md                  | 22 ++++++++++++++++++++++
 map2loop/version.py           |  2 +-
 3 files changed, 24 insertions(+), 2 deletions(-)

diff --git a/.release-please-manifest.json b/.release-please-manifest.json
index bbadce2a..c78fcfe9 100644
--- a/.release-please-manifest.json
+++ b/.release-please-manifest.json
@@ -1,3 +1,3 @@
 {
-  ".": "3.2.9"
+  ".": "3.2.10"
 }
\ No newline at end of file
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2729cc98..add53866 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,27 @@
 # Changelog
 
+## [3.2.10](https://github.com/Loop3D/map2loop/compare/v3.2.9...v3.2.10) (2026-01-20)
+
+
+### Bug Fixes
+
+* add ability to specify fault field to topoloy ([92f754b](https://github.com/Loop3D/map2loop/commit/92f754b2511a18f41c0f990ca6f6391774858ad1))
+* add call methods to allow for generic debug exporting ([3c0e0e0](https://github.com/Loop3D/map2loop/commit/3c0e0e075b9d91bc32341cc8e2f249a828b00ee4))
+* add option to set logger level and handler. ([5fb0720](https://github.com/Loop3D/map2loop/commit/5fb0720d795a7123094c39f171c3f497f53eda07))
+* adding safeguarding and fixing merge ([a675be2](https://github.com/Loop3D/map2loop/commit/a675be2d6cef5646a682dbcd5c706681830d1ad9))
+* applying copilot suggestions ([4021640](https://github.com/Loop3D/map2loop/commit/40216402634b5b9349719b681fb48525080b48a8))
+* beartype issues ([6e2ee69](https://github.com/Loop3D/map2loop/commit/6e2ee69f08d913ce3a55df3f08144ded01922829))
+* change to no default arguments and add unit_name_column as a required ([65a7b86](https://github.com/Loop3D/map2loop/commit/65a7b86ae94c44c098aceb17b13bf3fa624716b9))
+* only store geodataframe if layer is not empty ([46fade8](https://github.com/Loop3D/map2loop/commit/46fade8bc8a673c0b8c067115970da70fed209bc))
+* project import for test ([6eb1518](https://github.com/Loop3D/map2loop/commit/6eb151882369e3e7a062f60722aac9b7e3436736))
+* replace args for kwargs ([e1cb1e6](https://github.com/Loop3D/map2loop/commit/e1cb1e642ffc9a3c206d51a8ccf78a064d645b5f))
+* save location tracking as a line geodataframe for qgis compatibility ([17e14ca](https://github.com/Loop3D/map2loop/commit/17e14cae5cd466e8e3185e2569c311a611e0bed4))
+* set age based sorter args ([5226b49](https://github.com/Loop3D/map2loop/commit/5226b49cfeda5592d878fc07e5ff960ba125f51e))
+* update maximise contacts to compute length and check for empty contacts ([f2cf211](https://github.com/Loop3D/map2loop/commit/f2cf211ce61e9547ac9a00e4d70f00cba24f3d4d))
+* updating incorrect type hints ([71f473c](https://github.com/Loop3D/map2loop/commit/71f473cd7e18d984b474de3a1b4bfd91bac98d81))
+* updating logic ([86b658c](https://github.com/Loop3D/map2loop/commit/86b658c582c8ca6102729e8f8fd7370f28a113f3))
+* updating requirements for alpha sorter plus add check for empty contacts ([b0f12b6](https://github.com/Loop3D/map2loop/commit/b0f12b68cfea2e0d16935720f01240230d912595))
+
 ## [3.2.9](https://github.com/Loop3D/map2loop/compare/v3.2.8...v3.2.9) (2025-12-01)
 
 
diff --git a/map2loop/version.py b/map2loop/version.py
index b57bdcdb..149fe9ad 100644
--- a/map2loop/version.py
+++ b/map2loop/version.py
@@ -1 +1 @@
-__version__ = "3.2.9"
+__version__ = "3.2.10"

From 8cbb6fcb4f8d3a8b0d15c39f5e27f5b8e861fe61 Mon Sep 17 00:00:00 2001
From: rabii-chaarani <rabii.chaarani@monash.edu>
Date: Thu, 22 Jan 2026 10:30:16 +0930
Subject: [PATCH 29/31] fix: added utils imports

---
 map2loop/utils/__init__.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/map2loop/utils/__init__.py b/map2loop/utils/__init__.py
index c7fa49d6..42c4aa09 100644
--- a/map2loop/utils/__init__.py
+++ b/map2loop/utils/__init__.py
@@ -15,4 +15,8 @@
     normal_vector_to_dipdirection_dip,
     strike_dip_vector,
     generate_grid,
+    segment_measure_range,
+    clean_line_geometry,
+    nearest_orientation_to_line,
+    iter_line_segments
 )

From c4fe8809501c6e2d9d83fd862071b7b0e409308b Mon Sep 17 00:00:00 2001
From: rabii-chaarani <rabii.chaarani@monash.edu>
Date: Thu, 22 Jan 2026 10:29:50 +0930
Subject: [PATCH 30/31] fix: merged changes from master

---
 map2loop/utils/utility_functions.py | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/map2loop/utils/utility_functions.py b/map2loop/utils/utility_functions.py
index 42287d28..7e76ac64 100644
--- a/map2loop/utils/utility_functions.py
+++ b/map2loop/utils/utility_functions.py
@@ -591,7 +591,6 @@ def set_z_values_from_raster_df(dtm_data: gdal.Dataset, df: pandas.DataFrame):
     )
 
     return df
-<<<<<<< HEAD
 
 def clean_line_geometry(geometry: shapely.geometry.base.BaseGeometry):
     """Clean and normalize a Shapely geometry to a single LineString or None.
@@ -801,5 +800,3 @@ def segment_measure_range(parent_line: shapely.geometry.LineString, segment: sha
     if end_measure < start_measure:
         start_measure, end_measure = end_measure, start_measure
     return start_measure, end_measure
-=======
->>>>>>> 1febaa9 (style: formatting)

From 92ce07feebedd345bf53d74ec1cd03094a61486b Mon Sep 17 00:00:00 2001
From: rabii-chaarani <rabii.chaarani@monash.edu>
Date: Thu, 22 Jan 2026 12:21:44 +0930
Subject: [PATCH 31/31] fix: minor fixes

---
 map2loop/data_checks.py          | 2 +-
 map2loop/thickness_calculator.py | 1 -
 2 files changed, 1 insertion(+), 2 deletions(-)

diff --git a/map2loop/data_checks.py b/map2loop/data_checks.py
index 5d766267..5f031be0 100644
--- a/map2loop/data_checks.py
+++ b/map2loop/data_checks.py
@@ -460,7 +460,7 @@ def check_keys(d: dict, parent_key=""):
         "structure": {
             "orientation_type", "dipdir_column", "dip_column",
             "description_column", "bedding_text", "overturned_column", "overturned_text",
-            "objectid_column", "desciption_column",
+            "objectid_column", "desciption_column", "interp_source_column"
         },
         "geology": {
             "unitname_column", "alt_unitname_column", "group_column",
diff --git a/map2loop/thickness_calculator.py b/map2loop/thickness_calculator.py
index 5f6e3b57..a57632db 100644
--- a/map2loop/thickness_calculator.py
+++ b/map2loop/thickness_calculator.py
@@ -898,7 +898,6 @@ def compute(
             logger.error("Failed to create spatial index for geology data: %s", e, exc_info=True)
             geology_sindex = None
 
-        #TODO: check if sections and geology have same crs and reproject if needed. then check if sections overlap geology
         sections_bounds = sections.total_bounds
         geology_bounds = geology.total_bounds
         if (