Skip to content

Fix/interpolated thicknesses #187

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
wants to merge 6 commits into
base: master
Choose a base branch
from
Open

Fix/interpolated thicknesses #187

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
4e2e799
fix: accept pathlib.Path for save_path
rabii-chaarani Apr 23, 2025
ff23f45
fix: allow str and pathlib.Path
lachlangrose May 1, 2025
05bfb96
fix: interpolated thicknesses out by 1
lachlangrose May 1, 2025
66be5b0
Store using correct index
lachlangrose May 1, 2025
6a29f5d
fix: handle NaN values
rabii-chaarani Jun 5, 2025
899f699
fix: replace logging reference with logger
rabii-chaarani Jun 5, 2025
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
8 changes: 5 additions & 3 deletions map2loop/project.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1012,7 +1012,7 @@ def draw_geology_map(self, points: pandas.DataFrame = None, overlay: str = ""):
gdf.plot(ax=base, marker="o", color="red", markersize=5)

@beartype.beartype
def save_mapdata_to_files(self, save_path: str = ".", extension: str = ".shp.zip"):
def save_mapdata_to_files(self, save_path: Union[pathlib.Path,str], extension: str = ".shp.zip"):
"""
Saves the map data frames to csv files

Expand All @@ -1022,8 +1022,10 @@ def save_mapdata_to_files(self, save_path: str = ".", extension: str = ".shp.zip
extension (str, optional):
An alternate extension to save the GeoDataFrame in. Defaults to ".csv".
"""
if not os.path.exists(save_path):
os.mkdir(save_path)

save_path=pathlib.Path(save_path)
if not save_path.exists():
os.makedirs(save_path)
self.map_data.save_all_map_data(save_path, extension)

@beartype.beartype
Expand Down
63 changes: 38 additions & 25 deletions map2loop/thickness_calculator.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -317,16 +317,16 @@ def compute(
_dips = []
_location_tracking = []

for i in range(0, len(stratigraphic_order) - 1):
for i in (range(1, len(stratigraphic_order) )):
if (
stratigraphic_order[i] in basal_unit_list
and stratigraphic_order[i + 1] in basal_unit_list
and stratigraphic_order[i - 1] in basal_unit_list
):
basal_contact = contacts.loc[
contacts["basal_unit"] == stratigraphic_order[i]
contacts["basal_unit"] == stratigraphic_order[i-1]
].copy()
top_contact = basal_contacts.loc[
basal_contacts["basal_unit"] == stratigraphic_order[i + 1]
basal_contacts["basal_unit"] == stratigraphic_order[i]
].copy()
top_contact_geometry = [
shapely.geometry.shape(geom.__geo_interface__) for geom in top_contact.geometry
Expand Down Expand Up @@ -364,14 +364,18 @@ def compute(
p2[2] = map_data.get_value_from_raster(Datatype.DTM, p2[0], p2[1])
# calculate the length of the shortest line
line_length = scipy.spatial.distance.euclidean(p1, p2)
# find the indices of the points that are within 5% of the length of the shortest line
indices = shapely.dwithin(short_line, interp_points, line_length * 0.25)
# find the indices of the points that are within 10% of the length of the shortest line
indices = shapely.dwithin(short_line, interp_points, line_length * 1)
# get the dip of the points that are within
_dip = numpy.deg2rad(dip[indices])
_dips.append(_dip)
# calculate the true thickness t = L * sin(dip)
thickness = line_length * numpy.sin(_dip)

if all(numpy.isnan(dip[indices])):
pass
else:
_dip = numpy.nanmean(dip[indices])
_dip = numpy.deg2rad(_dip)
_dips.append(_dip)
# calculate the true thickness t = L * sin(dip)
thickness = line_length * numpy.sin(_dip)

# add location tracking
location_tracking = pandas.DataFrame(
{
Expand All @@ -384,28 +388,37 @@ def compute(
_location_tracking.append(location_tracking)

# Average thickness along the shortest line
if all(numpy.isnan(thickness)):
pass
else:
_thickness.append(numpy.nanmean(thickness))
_thickness.append(thickness)

# calculate the median thickness and standard deviation for the unit
_thickness = numpy.asarray(_thickness, dtype=numpy.float64)

median = numpy.nanmedian(_thickness)
mean = numpy.nanmean(_thickness)
std_dev = numpy.nanstd(_thickness, dtype=numpy.float64)

idx = thicknesses.index[
if all(numpy.isnan(_thickness)):
logger.warning(
f"Cannot calculate thickness of {stratigraphic_order[i + 1]}. Near dip data not found. Assign NaN to thickness."
)
idx = thicknesses.index[
thicknesses["name"] == stratigraphic_order[i + 1]
].tolist()[0]
thicknesses.loc[idx, "ThicknessMean"] = mean
thicknesses.loc[idx, "ThicknessMedian"] = median
thicknesses.loc[idx, "ThicknessStdDev"] = std_dev
thicknesses.loc[idx, "ThicknessMean"] = numpy.nan
thicknesses.loc[idx, "ThicknessMedian"] = numpy.nan
thicknesses.loc[idx, "ThicknessStdDev"] = numpy.nan

else:
median = numpy.nanmedian(_thickness)
mean = numpy.nanmean(_thickness)
std_dev = numpy.nanstd(_thickness, dtype=numpy.float64)

idx = thicknesses.index[
thicknesses["name"] == stratigraphic_order[i + 1]
].tolist()[0]
thicknesses.loc[idx, "ThicknessMean"] = mean
thicknesses.loc[idx, "ThicknessMedian"] = median
thicknesses.loc[idx, "ThicknessStdDev"] = std_dev

else:
logger.warning(
f"Thickness Calculator InterpolatedStructure: Cannot calculate thickness between {stratigraphic_order[i]} and {stratigraphic_order[i + 1]}\n"
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
Expand Down
Loading