Skip to content

DTM is optional #234

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

Merged
merged 3 commits into from
Sep 19, 2025
Merged

DTM is optional #234

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
510c6fc
fix: handle optional DTM data
rabii-chaarani Sep 17, 2025
46b9fb2
fix: update geometry after setting Z
rabii-chaarani Sep 17, 2025
01d08ed
fix: Z is calculated when DTM is available
rabii-chaarani Sep 17, 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
50 changes: 31 additions & 19 deletions map2loop/thickness_calculator.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -287,15 +287,24 @@ 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
contacts = set_z_values_from_raster_df(self.dtm_data, contacts)
# update the geometry of the contact points to include the Z value
contacts["geometry"] = contacts.apply(
lambda row: shapely.geometry.Point(row.geometry.x, row.geometry.y, row["Z"]), axis=1
)
contacts = set_z_values_from_raster_df(self.dtm_data, contacts)
# update the geometry of the contact points to include the Z value
contacts["geometry"] = contacts.apply(
lambda row: shapely.geometry.Point(row.geometry.x, row.geometry.y, row["Z"]), axis=1
)
else:
contacts["geometry"] = contacts.apply(
lambda row: shapely.geometry.Point(row.geometry.x, row.geometry.y,), axis=1
)
# spatial join the contact points with the basal contacts to get the unit for each contact point
contacts = contacts.sjoin(basal_contacts, how="inner", predicate="intersects")
contacts = contacts[["X", "Y", "Z", "geometry", "basal_unit"]].copy()
# keep only necessary columns
if 'Z' not in contacts.columns:
contacts = contacts[["X", "Y", "geometry", "basal_unit"]].copy()
if 'Z' in contacts.columns:
contacts = contacts[["X", "Y", "Z", "geometry", "basal_unit"]].copy()
# Interpolate the dip of the contacts
interpolator = DipDipDirectionInterpolator(data_type="dip")
# Interpolate the dip of the contacts
Expand All @@ -313,12 +322,13 @@ def compute(
interpolated_orientations.set_geometry(create_points(xy), inplace=True)
# set the crs of the interpolated orientations to the crs of the units
interpolated_orientations = interpolated_orientations.set_crs(crs=basal_contacts.crs)
# get the elevation Z of the interpolated points
interpolated = set_z_values_from_raster_df(self.dtm_data, interpolated_orientations)
# update the geometry of the interpolated points to include the Z value
interpolated["geometry"] = interpolated.apply(
lambda row: shapely.geometry.Point(row.geometry.x, row.geometry.y, row["Z"]), axis=1
)
if self.dtm_data is not None:
# get the elevation Z of the interpolated points
interpolated_orientations = set_z_values_from_raster_df(self.dtm_data, interpolated_orientations)
# update the geometry of the interpolated points to include the Z value
interpolated_orientations["geometry"] = interpolated_orientations.apply(
lambda row: shapely.geometry.Point(row.geometry.x, row.geometry.y, row["Z"]), axis=1
)
# for each interpolated point, assign name of unit using spatial join
units = geology_data.copy()
interpolated_orientations = interpolated_orientations.sjoin(
Expand Down Expand Up @@ -364,22 +374,24 @@ def compute(
# check if the short line is
if self.max_line_length is not None and short_line.length > self.max_line_length:
continue

inv_geotransform = gdal.InvGeoTransform(self.dtm_data.GetGeoTransform())
data_array = numpy.array(self.dtm_data.GetRasterBand(1).ReadAsArray().T)
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])
# get the elevation Z of the end point p1
p1[2] = value_from_raster(inv_geotransform, data_array, p1[0], p1[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])
# get the elevation Z of the end point p2
p2[2] = value_from_raster(inv_geotransform, data_array, p2[0], p2[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])
# 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
Expand Down
Loading