From 2cf8942c6e054e4dbb0a6ad4f5bc34f02216712c Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Mon, 21 Sep 2026 18:13:55 +0300 Subject: [PATCH] fix: correct InterpolatedStructure contact/dip selection and StructuralPoint strike wraparound InterpolatedStructure: - shapely.shortest_line(point, top_contact_geometry) broadcasts over the list of top-contact geometries and returns one line per entry. The code took short_line[0] unconditionally, so whenever a basal contact was split into more than one feature (common across faults or map-sheet edges) the thickness estimate used an arbitrary segment instead of the nearest one. Now the candidate with the smallest length is selected explicitly. - Dip was sampled from stratigraphic_order[i] (the overlying unit) while the measured line spans stratigraphic_order[i + 1] (the unit whose thickness is being computed). Dip is now sampled from the same unit being measured. StructuralPoint: - The strike-allowance check compared strike1/strike2 against a plain numeric range (strike - allowance, strike + allowance). Strike is a compass bearing that wraps at 360 degrees, so measurements near due north (e.g. strike=5, contact strike=350) were incorrectly rejected even though they are well within tolerance. Replaced with an angular difference check that handles the wraparound correctly. --- map2loop/thickness_calculator.py | 56 +++++++++++++++++++++++--------- 1 file changed, 41 insertions(+), 15 deletions(-) diff --git a/map2loop/thickness_calculator.py b/map2loop/thickness_calculator.py index 483975d0..8280fcb8 100644 --- a/map2loop/thickness_calculator.py +++ b/map2loop/thickness_calculator.py @@ -34,6 +34,16 @@ from osgeo import gdal from shapely.errors import UnsupportedGEOSVersionError + +def _angular_difference(angle_a: float, angle_b: float) -> float: + """ + Smallest difference between two compass bearings (0-360 degrees), correctly + handling wraparound e.g. the difference between 350 and 5 degrees is 15, not 345. + """ + diff = abs(angle_a - angle_b) % 360 + return min(diff, 360 - diff) + + class ThicknessCalculator(ABC): """ Base Class of Thickness Calculator used to force structure of ThicknessCalculator @@ -368,20 +378,32 @@ def compute( shapely.geometry.shape(geom.__geo_interface__) for geom in top_contact.geometry ] if basal_contact is not None and top_contact is not None: + # dip is sampled from the unit whose thickness is being measured + # (stratigraphic_order[i + 1], bounded above by basal_contact and + # below by top_contact), not from the overlying unit. interp_points = interpolated_orientations.loc[ - interpolated_orientations["UNITNAME"] == stratigraphic_order[i], "geometry" + interpolated_orientations["UNITNAME"] == stratigraphic_order[i + 1], "geometry" ].copy() dip = interpolated_orientations.loc[ - interpolated_orientations["UNITNAME"] == stratigraphic_order[i], "dip" + interpolated_orientations["UNITNAME"] == stratigraphic_order[i + 1], "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) + # find the shortest line between the basal contact point and every + # top contact geometry, then keep the globally shortest one. shapely + # broadcasts a scalar point against a list of geometries and returns + # one line per list entry, so a single geometry cannot simply be + # indexed out without checking which candidate is actually nearest. + short_line_candidates = numpy.atleast_1d( + shapely.shortest_line(row.geometry, top_contact_geometry) + ) + short_line = short_line_candidates[ + numpy.argmin(shapely.length(short_line_candidates)) + ] # check if the short line is - if self.max_line_length is not None and short_line[0].length > self.max_line_length: + 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()) @@ -389,15 +411,15 @@ def compute( # extract the end points of the shortest line p1 = numpy.zeros(3) - p1[0] = numpy.asarray(short_line[0].coords[0][0]) - p1[1] = numpy.asarray(short_line[0].coords[0][1]) + p1[0] = numpy.asarray(short_line.coords[0][0]) + p1[1] = numpy.asarray(short_line.coords[0][1]) if self.dtm_data is not None: # get the elevation Z of the end point p1 p1[2] = value_from_raster(inv_geotransform, data_array, p1[0], p1[1]) # create array to store xyz coordinates of the end point p2 p2 = numpy.zeros(3) - p2[0] = numpy.asarray(short_line[0].coords[-1][0]) - p2[1] = numpy.asarray(short_line[0].coords[-1][1]) + p2[0] = numpy.asarray(short_line.coords[-1][0]) + p2[1] = numpy.asarray(short_line.coords[-1][1]) if self.dtm_data is not None: # get the elevation Z of the end point p2 p2[2] = value_from_raster(inv_geotransform, data_array, p2[0], p2[1]) @@ -406,13 +428,13 @@ def compute( # find the indices of the points that are within 5% of the length of the shortest line try: # GEOS 3.10.0+ - indices = shapely.dwithin(short_line[0], interp_points, line_length * 0.25) + indices = shapely.dwithin(short_line, interp_points, line_length * 0.25) except UnsupportedGEOSVersionError: - indices= numpy.array([shapely.distance(short_line[0],point)<= (line_length * 0.25) for point in interp_points]) + indices= numpy.array([shapely.distance(short_line,point)<= (line_length * 0.25) for point in interp_points]) # get the dip of the points that are within _dip = numpy.deg2rad(dip[indices]) if len(_dip) > 0: - _lines.extend([short_line[0]]*len(_dip)) + _lines.extend([short_line]*len(_dip)) _dips.extend(_dip) # calculate the true thickness t = L * sin(dip) thickness = line_length * numpy.sin(_dip) @@ -699,9 +721,13 @@ def compute( strike1 = find_segment_strike_from_pt(seg1, int_pt1, measurement) strike2 = find_segment_strike_from_pt(seg2, int_pt2, measurement) - # check to see if the strike of the stratigraphic measurement is within the strike allowance of the strike of the geological contact - b_s = strike - self.strike_allowance, strike + self.strike_allowance - if not (b_s[0] < strike1 < b_s[1] and b_s[0] < strike2 < b_s[1]): + # check to see if the strike of the stratigraphic measurement is within the strike allowance + # of the strike of the geological contact. Strike is a compass bearing (wraps at 360 degrees), + # so the comparison must use angular difference rather than a plain numeric range. + if ( + _angular_difference(strike, strike1) > self.strike_allowance + or _angular_difference(strike, strike2) > self.strike_allowance + ): continue # build the debug info