Skip to content

Commit d5e5862

Browse files
noellehmchengnoellehmcheng
authored
Merge pull request #25 from Loop3D/noelle/thickness_calculator
thickness calculator
2 parents d4d6d83 + 291de7a commit d5e5862

File tree

2 files changed

+124
-36
lines changed

2 files changed

+124
-36
lines changed

‎m2l/processing/algorithms/sampler.py‎

Lines changed: 28 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -176,11 +176,19 @@ def processAlgorithm(
176176
samples = sampler.sample(spatial_data_gdf)
177177

178178
fields = QgsFields()
179-
fields.append(QgsField("ID", QVariant.String))
180-
fields.append(QgsField("X", QVariant.Double))
181-
fields.append(QgsField("Y", QVariant.Double))
182-
fields.append(QgsField("Z", QVariant.Double))
183-
fields.append(QgsField("featureId", QVariant.String))
179+
if samples is not None and not samples.empty:
180+
for column_name in samples.columns:
181+
dtype = samples[column_name].dtype
182+
dtype_str = str(dtype)
183+
184+
if dtype_str in ['float16', 'float32', 'float64']:
185+
field_type = QVariant.Double
186+
elif dtype_str in ['int8', 'int16', 'int32', 'int64']:
187+
field_type = QVariant.Int
188+
else:
189+
field_type = QVariant.String
190+
191+
fields.append(QgsField(column_name, field_type))
184192

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

210-
feature.setAttributes([
211-
str(row.get('ID', '')),
212-
float(row.get('X', 0)),
213-
float(row.get('Y', 0)),
214-
float(row.get('Z', 0)) if pd.notna(row.get('Z')) else 0.0,
215-
str(row.get('featureId', ''))
216-
])
218+
attributes = []
219+
for column_name in samples.columns:
220+
value = row.get(column_name)
221+
dtype = samples[column_name].dtype
222+
223+
if pd.isna(value):
224+
attributes.append(None)
225+
elif dtype in ['float16', 'float32', 'float64']:
226+
attributes.append(float(value))
227+
elif dtype in ['int8', 'int16', 'int32', 'int64']:
228+
attributes.append(int(value))
229+
else:
230+
attributes.append(str(value))
231+
232+
feature.setAttributes(attributes)
217233

218234
sink.addFeature(feature)
219235

‎m2l/processing/algorithms/thickness_calculator.py‎

Lines changed: 96 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@
1010
"""
1111
# Python imports
1212
from typing import Any, Optional
13+
import pandas as pd
1314

1415
# QGIS imports
1516
from qgis import processing
@@ -28,7 +29,6 @@
2829
QgsProcessingParameterMatrix,
2930
QgsSettings,
3031
QgsProcessingParameterRasterLayer,
31-
QgsProcessingParameterMapLayer
3232
)
3333
# Internal imports
3434
from ...main.vectorLayerWrapper import (
@@ -48,6 +48,7 @@ class ThicknessCalculatorAlgorithm(QgsProcessingAlgorithm):
4848

4949
INPUT_THICKNESS_CALCULATOR_TYPE = 'THICKNESS_CALCULATOR_TYPE'
5050
INPUT_DTM = 'DTM'
51+
INPUT_BOUNDING_BOX_TYPE = 'BOUNDING_BOX_TYPE'
5152
INPUT_BOUNDING_BOX = 'BOUNDING_BOX'
5253
INPUT_MAX_LINE_LENGTH = 'MAX_LINE_LENGTH'
5354
INPUT_STRATI_COLUMN = 'STRATIGRAPHIC_COLUMN'
@@ -56,8 +57,10 @@ class ThicknessCalculatorAlgorithm(QgsProcessingAlgorithm):
5657
INPUT_DIPDIR_FIELD = 'DIPDIR_FIELD'
5758
INPUT_DIP_FIELD = 'DIP_FIELD'
5859
INPUT_GEOLOGY = 'GEOLOGY'
60+
INPUT_ORIENTATION_TYPE = 'ORIENTATION_TYPE'
5961
INPUT_UNIT_NAME_FIELD = 'UNIT_NAME_FIELD'
6062
INPUT_SAMPLED_CONTACTS = 'SAMPLED_CONTACTS'
63+
INPUT_STRATIGRAPHIC_COLUMN_LAYER = 'STRATIGRAPHIC_COLUMN_LAYER'
6164

6265
OUTPUT = "THICKNESS"
6366

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

104+
self.addParameter(
105+
QgsProcessingParameterEnum(
106+
self.INPUT_BOUNDING_BOX_TYPE,
107+
"Bounding Box Type",
108+
options=['Extract from geology layer', 'User defined'],
109+
allowMultiple=False,
110+
defaultValue=1
111+
)
112+
)
113+
101114
bbox_settings = QgsSettings()
102115
last_bbox = bbox_settings.value("m2l/bounding_box", "")
103116
self.addParameter(
104117
QgsProcessingParameterMatrix(
105118
self.INPUT_BOUNDING_BOX,
106-
description="Bounding Box",
119+
description="Static Bounding Box",
107120
headers=['minx','miny','maxx','maxy'],
108121
numberRows=1,
109-
defaultValue=last_bbox
122+
defaultValue=last_bbox,
123+
optional=True
110124
)
111125
)
126+
112127
self.addParameter(
113128
QgsProcessingParameterNumber(
114129
self.INPUT_MAX_LINE_LENGTH,
@@ -142,13 +157,26 @@ def initAlgorithm(self, config: Optional[dict[str, Any]] = None) -> None:
142157
defaultValue='Formation'
143158
)
144159
)
145-
160+
146161
self.addParameter(
147162
QgsProcessingParameterFeatureSource(
148-
self.INPUT_STRATI_COLUMN,
149-
"Stratigraphic Order",
163+
'STRATIGRAPHIC_COLUMN_LAYER',
164+
'Stratigraphic Column Layer (from sorter)',
150165
[QgsProcessing.TypeVector],
151-
defaultValue='formation',
166+
optional=True
167+
)
168+
)
169+
170+
strati_settings = QgsSettings()
171+
last_strati_column = strati_settings.value("m2l/strati_column", "")
172+
self.addParameter(
173+
QgsProcessingParameterMatrix(
174+
name=self.INPUT_STRATI_COLUMN,
175+
description="Stratigraphic Order",
176+
headers=["Unit"],
177+
numberRows=0,
178+
defaultValue=last_strati_column,
179+
optional=True
152180
)
153181
)
154182
self.addParameter(
@@ -165,6 +193,14 @@ def initAlgorithm(self, config: Optional[dict[str, Any]] = None) -> None:
165193
[QgsProcessing.TypeVectorPoint],
166194
)
167195
)
196+
self.addParameter(
197+
QgsProcessingParameterEnum(
198+
self.INPUT_ORIENTATION_TYPE,
199+
'Orientation Type',
200+
options=['Dip Direction', 'Strike'],
201+
defaultValue=0 # Default to Dip Direction
202+
)
203+
)
168204
self.addParameter(
169205
QgsProcessingParameterField(
170206
self.INPUT_DIPDIR_FIELD,
@@ -201,27 +237,60 @@ def processAlgorithm(
201237
thickness_type_index = self.parameterAsEnum(parameters, self.INPUT_THICKNESS_CALCULATOR_TYPE, context)
202238
thickness_type = ['InterpolatedStructure', 'StructuralPoint'][thickness_type_index]
203239
dtm_data = self.parameterAsRasterLayer(parameters, self.INPUT_DTM, context)
204-
bounding_box = self.parameterAsMatrix(parameters, self.INPUT_BOUNDING_BOX, context)
240+
bounding_box_type = self.parameterAsEnum(parameters, self.INPUT_BOUNDING_BOX_TYPE, context)
205241
max_line_length = self.parameterAsSource(parameters, self.INPUT_MAX_LINE_LENGTH, context)
206242
basal_contacts = self.parameterAsSource(parameters, self.INPUT_BASAL_CONTACTS, context)
207243
geology_data = self.parameterAsSource(parameters, self.INPUT_GEOLOGY, context)
208-
stratigraphic_order = self.parameterAsSource(parameters, self.INPUT_STRATI_COLUMN, context)
209244
structure_data = self.parameterAsSource(parameters, self.INPUT_STRUCTURE_DATA, context)
245+
orientation_type = self.parameterAsEnum(parameters, self.INPUT_ORIENTATION_TYPE, context)
246+
is_strike = (orientation_type == 1)
210247
structure_dipdir_field = self.parameterAsString(parameters, self.INPUT_DIPDIR_FIELD, context)
211248
structure_dip_field = self.parameterAsString(parameters, self.INPUT_DIP_FIELD, context)
212249
sampled_contacts = self.parameterAsSource(parameters, self.INPUT_SAMPLED_CONTACTS, context)
213250
unit_name_field = self.parameterAsString(parameters, self.INPUT_UNIT_NAME_FIELD, context)
214251

215-
bbox_settings = QgsSettings()
216-
bbox_settings.setValue("m2l/bounding_box", bounding_box)
217-
218-
if isinstance(stratigraphic_order, QgsProcessingParameterMapLayer) :
219-
raise QgsProcessingException("Invalid stratigraphic column layer")
220-
221-
if stratigraphic_order is not None:
222-
# extract unit names from stratigraphic_order
223-
field_name = "unit_name"
224-
stratigraphic_order = [f[field_name] for f in stratigraphic_order.getFeatures()]
252+
if bounding_box_type == 0:
253+
geology_layer = self.parameterAsVectorLayer(parameters, self.INPUT_GEOLOGY, context)
254+
extent = geology_layer.extent()
255+
bounding_box = {
256+
'minx': extent.xMinimum(),
257+
'miny': extent.yMinimum(),
258+
'maxx': extent.xMaximum(),
259+
'maxy': extent.yMaximum()
260+
}
261+
feedback.pushInfo("Using bounding box from geology layer")
262+
else:
263+
static_bbox_matrix = self.parameterAsMatrix(parameters, self.INPUT_BOUNDING_BOX, context)
264+
if not static_bbox_matrix or len(static_bbox_matrix) == 0:
265+
raise QgsProcessingException("Bounding box is required")
266+
267+
bounding_box = matrixToDict(static_bbox_matrix)
268+
269+
bbox_settings = QgsSettings()
270+
bbox_settings.setValue("m2l/bounding_box", static_bbox_matrix)
271+
feedback.pushInfo("Using bounding box from user input")
272+
273+
stratigraphic_column_source = self.parameterAsSource(parameters, self.INPUT_STRATIGRAPHIC_COLUMN_LAYER, context)
274+
stratigraphic_order = []
275+
if stratigraphic_column_source is not None:
276+
ordered_pairs =[]
277+
for feature in stratigraphic_column_source.getFeatures():
278+
order = feature.attribute('order')
279+
unit_name = feature.attribute('unit_name')
280+
ordered_pairs.append((order, unit_name))
281+
ordered_pairs.sort(key=lambda x: x[0])
282+
stratigraphic_order = [pair[1] for pair in ordered_pairs]
283+
feedback.pushInfo(f"DEBUG: parameterAsVectorLayer Stratigraphic order: {stratigraphic_order}")
284+
else:
285+
matrix_stratigraphic_order = self.parameterAsMatrix(parameters, self.INPUT_STRATI_COLUMN, context)
286+
if matrix_stratigraphic_order:
287+
stratigraphic_order = [str(row) for row in matrix_stratigraphic_order if row]
288+
else:
289+
raise QgsProcessingException("Stratigraphic column layer is required")
290+
if stratigraphic_order:
291+
matrix = [[unit] for unit in stratigraphic_order]
292+
strati_column_settings = QgsSettings()
293+
strati_column_settings.setValue('m2l/strati_column', matrix)
225294
# convert layers to dataframe or geodataframe
226295
units = qgsLayerToDataFrame(geology_data)
227296
geology_data = qgsLayerToGeoDataFrame(geology_data)
@@ -231,9 +300,8 @@ def processAlgorithm(
231300
missing_fields = []
232301
if unit_name_field != 'UNITNAME' and unit_name_field in geology_data.columns:
233302
geology_data = geology_data.rename(columns={unit_name_field: 'UNITNAME'})
234-
units = units.rename(columns={unit_name_field: 'UNITNAME'})
235-
units = units.drop_duplicates(subset=['UNITNAME']).reset_index(drop=True)
236-
units = units.rename(columns={'UNITNAME': 'name'})
303+
units_unique = units.drop_duplicates(subset=[unit_name_field]).reset_index(drop=True)
304+
units = pd.DataFrame({'name': units_unique[unit_name_field]})
237305
if structure_data is not None:
238306
if structure_dipdir_field:
239307
if structure_dipdir_field in structure_data.columns:
@@ -251,14 +319,17 @@ def processAlgorithm(
251319
)
252320
if rename_map:
253321
structure_data = structure_data.rename(columns=rename_map)
322+
254323
sampled_contacts = qgsLayerToDataFrame(sampled_contacts)
324+
sampled_contacts['X'] = sampled_contacts['X'].astype(float)
325+
sampled_contacts['Y'] = sampled_contacts['Y'].astype(float)
326+
sampled_contacts['Z'] = sampled_contacts['Z'].astype(float)
255327
dtm_data = qgsRasterToGdalDataset(dtm_data)
256-
bounding_box = matrixToDict(bounding_box)
257-
feedback.pushInfo("Calculating unit thicknesses...")
258328
if thickness_type == "InterpolatedStructure":
259329
thickness_calculator = InterpolatedStructure(
260330
dtm_data=dtm_data,
261331
bounding_box=bounding_box,
332+
is_strike=is_strike
262333
)
263334
thicknesses = thickness_calculator.compute(
264335
units,
@@ -274,6 +345,7 @@ def processAlgorithm(
274345
dtm_data=dtm_data,
275346
bounding_box=bounding_box,
276347
max_line_length=max_line_length,
348+
is_strike=is_strike
277349
)
278350
thicknesses =thickness_calculator.compute(
279351
units,

0 commit comments

Comments
 (0)