Skip to content

thickness calculator #25

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
noellehmcheng merged 20 commits into from
Oct 16, 2025
Merged

thickness calculator #25

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
20 commits
Select commit Hold shift + click to select a range
73b6e22
feat: thickness calculator tool
rabii-chaarani Sep 17, 2025
b483f00
feat: raster and dataframe handling
rabii-chaarani Sep 17, 2025
09091a4
Merge branch 'processing/processing_tools' into processing/thickness_…
rabii-chaarani Sep 22, 2025
90663cf
fix: remove unused QgsSettings
rabii-chaarani Sep 22, 2025
4ff7546
fix: updated tearDownClass
rabii-chaarani Sep 22, 2025
14db8d0
fix: add pass statement
rabii-chaarani Sep 22, 2025
ec993bb
update sampler tests
noellehmcheng Sep 22, 2025
83c1703
update sampler
noellehmcheng Sep 22, 2025
747abef
Merge branch 'processing/thickness_calculator' of https://github.com/…
noellehmcheng Sep 22, 2025
5205cf5
update strati column and bounding box
noellehmcheng Sep 23, 2025
a259cae
fix sample contacts data type
noellehmcheng Sep 23, 2025
00d4bb8
add orientation type for structure data
noellehmcheng Sep 23, 2025
f59f4f8
fix units in ThicknessCalculatorAlgorithm
noellehmcheng Sep 24, 2025
2cd4c53
feat dynamic field handling in Sampler
noellehmcheng Oct 6, 2025
97c5274
fix strat column data in ThicknessCalculator
noellehmcheng Oct 6, 2025
fe429df
merge with processing/processing_tools
noellehmcheng Oct 7, 2025
53917e3
add user defined boundingbox in ThicknessCalculatorAlgorithm
noellehmcheng Oct 15, 2025
c94b40c
replace QMetaType with QVariant in sampler
noellehmcheng Oct 15, 2025
b32baec
change orientation type variable name in ThicknessCalculatorAlgorithm
noellehmcheng Oct 15, 2025
291de7a
Merge branch 'processing/processing_tools' into noelle/thickness_calc…
noellehmcheng Oct 16, 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
40 changes: 28 additions & 12 deletions m2l/processing/algorithms/sampler.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -176,11 +176,19 @@ def processAlgorithm(
samples = sampler.sample(spatial_data_gdf)

fields = QgsFields()
fields.append(QgsField("ID", QVariant.String))
fields.append(QgsField("X", QVariant.Double))
Copy link
Member

@rabii-chaarani rabii-chaarani Oct 8, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should use QVariant for the moment, as QMetaType don't work for older versions

fields.append(QgsField("Y", QVariant.Double))
fields.append(QgsField("Z", QVariant.Double))
fields.append(QgsField("featureId", QVariant.String))
if samples is not None and not samples.empty:
for column_name in samples.columns:
dtype = samples[column_name].dtype
dtype_str = str(dtype)

if dtype_str in ['float16', 'float32', 'float64']:
field_type = QVariant.Double
elif dtype_str in ['int8', 'int16', 'int32', 'int64']:
field_type = QVariant.Int
else:
field_type = QVariant.String

fields.append(QgsField(column_name, field_type))

crs = None
if spatial_data_gdf is not None and spatial_data_gdf.crs is not None:
Expand All @@ -207,13 +215,21 @@ def processAlgorithm(
#spacing has no z values
feature.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(row['X'], row['Y'])))

feature.setAttributes([
str(row.get('ID', '')),
float(row.get('X', 0)),
float(row.get('Y', 0)),
float(row.get('Z', 0)) if pd.notna(row.get('Z')) else 0.0,
str(row.get('featureId', ''))
])
attributes = []
for column_name in samples.columns:
value = row.get(column_name)
dtype = samples[column_name].dtype

if pd.isna(value):
attributes.append(None)
elif dtype in ['float16', 'float32', 'float64']:
attributes.append(float(value))
elif dtype in ['int8', 'int16', 'int32', 'int64']:
attributes.append(int(value))
else:
attributes.append(str(value))

feature.setAttributes(attributes)

sink.addFeature(feature)

Expand Down
120 changes: 96 additions & 24 deletions m2l/processing/algorithms/thickness_calculator.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
"""
# Python imports
from typing import Any, Optional
import pandas as pd

# QGIS imports
from qgis import processing
Expand All @@ -28,7 +29,6 @@
QgsProcessingParameterMatrix,
QgsSettings,
QgsProcessingParameterRasterLayer,
QgsProcessingParameterMapLayer
)
# Internal imports
from ...main.vectorLayerWrapper import (
Expand All @@ -48,6 +48,7 @@ class ThicknessCalculatorAlgorithm(QgsProcessingAlgorithm):

INPUT_THICKNESS_CALCULATOR_TYPE = 'THICKNESS_CALCULATOR_TYPE'
INPUT_DTM = 'DTM'
INPUT_BOUNDING_BOX_TYPE = 'BOUNDING_BOX_TYPE'
INPUT_BOUNDING_BOX = 'BOUNDING_BOX'
INPUT_MAX_LINE_LENGTH = 'MAX_LINE_LENGTH'
INPUT_STRATI_COLUMN = 'STRATIGRAPHIC_COLUMN'
Expand All @@ -56,8 +57,10 @@ class ThicknessCalculatorAlgorithm(QgsProcessingAlgorithm):
INPUT_DIPDIR_FIELD = 'DIPDIR_FIELD'
INPUT_DIP_FIELD = 'DIP_FIELD'
INPUT_GEOLOGY = 'GEOLOGY'
INPUT_ORIENTATION_TYPE = 'ORIENTATION_TYPE'
INPUT_UNIT_NAME_FIELD = 'UNIT_NAME_FIELD'
INPUT_SAMPLED_CONTACTS = 'SAMPLED_CONTACTS'
INPUT_STRATIGRAPHIC_COLUMN_LAYER = 'STRATIGRAPHIC_COLUMN_LAYER'

OUTPUT = "THICKNESS"

Expand Down Expand Up @@ -98,17 +101,29 @@ def initAlgorithm(self, config: Optional[dict[str, Any]] = None) -> None:
)
)

self.addParameter(
QgsProcessingParameterEnum(
self.INPUT_BOUNDING_BOX_TYPE,
"Bounding Box Type",
options=['Extract from geology layer', 'User defined'],
allowMultiple=False,
defaultValue=1
)
)

bbox_settings = QgsSettings()
last_bbox = bbox_settings.value("m2l/bounding_box", "")
self.addParameter(
QgsProcessingParameterMatrix(
self.INPUT_BOUNDING_BOX,
description="Bounding Box",
description="Static Bounding Box",
headers=['minx','miny','maxx','maxy'],
numberRows=1,
defaultValue=last_bbox
defaultValue=last_bbox,
optional=True
)
)

self.addParameter(
QgsProcessingParameterNumber(
self.INPUT_MAX_LINE_LENGTH,
Expand Down Expand Up @@ -142,13 +157,26 @@ def initAlgorithm(self, config: Optional[dict[str, Any]] = None) -> None:
defaultValue='Formation'
)
)

self.addParameter(
QgsProcessingParameterFeatureSource(
self.INPUT_STRATI_COLUMN,
"Stratigraphic Order",
'STRATIGRAPHIC_COLUMN_LAYER',
'Stratigraphic Column Layer (from sorter)',
[QgsProcessing.TypeVector],
defaultValue='formation',
optional=True
)
)

strati_settings = QgsSettings()
last_strati_column = strati_settings.value("m2l/strati_column", "")
self.addParameter(
QgsProcessingParameterMatrix(
name=self.INPUT_STRATI_COLUMN,
description="Stratigraphic Order",
headers=["Unit"],
numberRows=0,
defaultValue=last_strati_column,
optional=True
)
)
self.addParameter(
Expand All @@ -165,6 +193,14 @@ def initAlgorithm(self, config: Optional[dict[str, Any]] = None) -> None:
[QgsProcessing.TypeVectorPoint],
)
)
self.addParameter(
QgsProcessingParameterEnum(
self.INPUT_ORIENTATION_TYPE,
'Orientation Type',
options=['Dip Direction', 'Strike'],
defaultValue=0 # Default to Dip Direction
)
)
self.addParameter(
QgsProcessingParameterField(
self.INPUT_DIPDIR_FIELD,
Expand Down Expand Up @@ -201,27 +237,60 @@ def processAlgorithm(
thickness_type_index = self.parameterAsEnum(parameters, self.INPUT_THICKNESS_CALCULATOR_TYPE, context)
thickness_type = ['InterpolatedStructure', 'StructuralPoint'][thickness_type_index]
dtm_data = self.parameterAsRasterLayer(parameters, self.INPUT_DTM, context)
bounding_box = self.parameterAsMatrix(parameters, self.INPUT_BOUNDING_BOX, context)
bounding_box_type = self.parameterAsEnum(parameters, self.INPUT_BOUNDING_BOX_TYPE, context)
max_line_length = self.parameterAsSource(parameters, self.INPUT_MAX_LINE_LENGTH, context)
basal_contacts = self.parameterAsSource(parameters, self.INPUT_BASAL_CONTACTS, context)
geology_data = self.parameterAsSource(parameters, self.INPUT_GEOLOGY, context)
stratigraphic_order = self.parameterAsSource(parameters, self.INPUT_STRATI_COLUMN, context)
structure_data = self.parameterAsSource(parameters, self.INPUT_STRUCTURE_DATA, context)
orientation_type = self.parameterAsEnum(parameters, self.INPUT_ORIENTATION_TYPE, context)
is_strike = (orientation_type == 1)
structure_dipdir_field = self.parameterAsString(parameters, self.INPUT_DIPDIR_FIELD, context)
structure_dip_field = self.parameterAsString(parameters, self.INPUT_DIP_FIELD, context)
sampled_contacts = self.parameterAsSource(parameters, self.INPUT_SAMPLED_CONTACTS, context)
unit_name_field = self.parameterAsString(parameters, self.INPUT_UNIT_NAME_FIELD, context)

bbox_settings = QgsSettings()
bbox_settings.setValue("m2l/bounding_box", bounding_box)

if isinstance(stratigraphic_order, QgsProcessingParameterMapLayer) :
raise QgsProcessingException("Invalid stratigraphic column layer")

if stratigraphic_order is not None:
# extract unit names from stratigraphic_order
field_name = "unit_name"
stratigraphic_order = [f[field_name] for f in stratigraphic_order.getFeatures()]
if bounding_box_type == 0:
geology_layer = self.parameterAsVectorLayer(parameters, self.INPUT_GEOLOGY, context)
extent = geology_layer.extent()
bounding_box = {
'minx': extent.xMinimum(),
'miny': extent.yMinimum(),
'maxx': extent.xMaximum(),
'maxy': extent.yMaximum()
}
feedback.pushInfo("Using bounding box from geology layer")
else:
static_bbox_matrix = self.parameterAsMatrix(parameters, self.INPUT_BOUNDING_BOX, context)
if not static_bbox_matrix or len(static_bbox_matrix) == 0:
raise QgsProcessingException("Bounding box is required")

bounding_box = matrixToDict(static_bbox_matrix)

bbox_settings = QgsSettings()
bbox_settings.setValue("m2l/bounding_box", static_bbox_matrix)
feedback.pushInfo("Using bounding box from user input")

stratigraphic_column_source = self.parameterAsSource(parameters, self.INPUT_STRATIGRAPHIC_COLUMN_LAYER, context)
stratigraphic_order = []
if stratigraphic_column_source is not None:
ordered_pairs =[]
for feature in stratigraphic_column_source.getFeatures():
order = feature.attribute('order')
unit_name = feature.attribute('unit_name')
ordered_pairs.append((order, unit_name))
ordered_pairs.sort(key=lambda x: x[0])
stratigraphic_order = [pair[1] for pair in ordered_pairs]
feedback.pushInfo(f"DEBUG: parameterAsVectorLayer Stratigraphic order: {stratigraphic_order}")
else:
matrix_stratigraphic_order = self.parameterAsMatrix(parameters, self.INPUT_STRATI_COLUMN, context)
if matrix_stratigraphic_order:
stratigraphic_order = [str(row) for row in matrix_stratigraphic_order if row]
else:
raise QgsProcessingException("Stratigraphic column layer is required")
if stratigraphic_order:
matrix = [[unit] for unit in stratigraphic_order]
strati_column_settings = QgsSettings()
strati_column_settings.setValue('m2l/strati_column', matrix)
# convert layers to dataframe or geodataframe
units = qgsLayerToDataFrame(geology_data)
geology_data = qgsLayerToGeoDataFrame(geology_data)
Expand All @@ -231,9 +300,8 @@ def processAlgorithm(
missing_fields = []
if unit_name_field != 'UNITNAME' and unit_name_field in geology_data.columns:
geology_data = geology_data.rename(columns={unit_name_field: 'UNITNAME'})
units = units.rename(columns={unit_name_field: 'UNITNAME'})
units = units.drop_duplicates(subset=['UNITNAME']).reset_index(drop=True)
units = units.rename(columns={'UNITNAME': 'name'})
units_unique = units.drop_duplicates(subset=[unit_name_field]).reset_index(drop=True)
units = pd.DataFrame({'name': units_unique[unit_name_field]})
if structure_data is not None:
if structure_dipdir_field:
if structure_dipdir_field in structure_data.columns:
Expand All @@ -251,14 +319,17 @@ def processAlgorithm(
)
if rename_map:
structure_data = structure_data.rename(columns=rename_map)

sampled_contacts = qgsLayerToDataFrame(sampled_contacts)
sampled_contacts['X'] = sampled_contacts['X'].astype(float)
sampled_contacts['Y'] = sampled_contacts['Y'].astype(float)
sampled_contacts['Z'] = sampled_contacts['Z'].astype(float)
dtm_data = qgsRasterToGdalDataset(dtm_data)
bounding_box = matrixToDict(bounding_box)
feedback.pushInfo("Calculating unit thicknesses...")
if thickness_type == "InterpolatedStructure":
thickness_calculator = InterpolatedStructure(
dtm_data=dtm_data,
bounding_box=bounding_box,
is_strike=is_strike
)
thicknesses = thickness_calculator.compute(
units,
Expand All @@ -274,6 +345,7 @@ def processAlgorithm(
dtm_data=dtm_data,
bounding_box=bounding_box,
max_line_length=max_line_length,
is_strike=is_strike
)
thicknesses =thickness_calculator.compute(
units,
Expand Down