Skip to content

fix(map2loop): InterpolatedStructure contact/dip selection and StructuralPoint strike wraparound #327

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
lachlangrose wants to merge 1 commit into
base: master
Choose a base branch
from
Open

fix(map2loop): InterpolatedStructure contact/dip selection and StructuralPoint strike wraparound #327

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
56 changes: 41 additions & 15 deletions packages/map2loop/src/map2loop/thickness_calculator.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -368,36 +378,48 @@ 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())
data_array = numpy.array(self.dtm_data.GetRasterBand(1).ReadAsArray().T)

# extract the end points of the shortest line
p1 = numpy.zeros(3)
p1[0] = numpy.asarray(short_line[0].coords[0][0])
p1[1] = numpy.asarray(short_line[0].coords[0][1])
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])
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand Down
Loading