|
12 | 12 | from map2loop.thickness_calculator import InterpolatedStructure, StructuralPoint |
13 | 13 | from osgeo import gdal |
14 | 14 |
|
15 | | -from ..main.vectorLayerWrapper import qgsLayerToGeoDataFrame |
| 15 | +from ..main.vectorLayerWrapper import qgsLayerToDataFrame, qgsLayerToGeoDataFrame |
16 | 16 |
|
17 | 17 | # Mapping of sorter names to sorter classes |
18 | 18 | SORTER_LIST = { |
@@ -340,8 +340,13 @@ def calculate_thickness( |
340 | 340 | geology_gdf = qgsLayerToGeoDataFrame(geology) |
341 | 341 | basal_contacts_gdf = qgsLayerToGeoDataFrame(basal_contacts) |
342 | 342 | sampled_contacts_gdf = qgsLayerToGeoDataFrame(sampled_contacts) |
343 | | - structure_gdf = qgsLayerToGeoDataFrame(structure) |
344 | | - |
| 343 | + structure_gdf = qgsLayerToDataFrame(structure) |
| 344 | + bounding_box = { |
| 345 | + 'maxx': geology_gdf.total_bounds[2], |
| 346 | + 'minx': geology_gdf.total_bounds[0], |
| 347 | + 'maxy': geology_gdf.total_bounds[3], |
| 348 | + 'miny': geology_gdf.total_bounds[1], |
| 349 | + } |
345 | 350 | # Rename unit name field if needed |
346 | 351 | if unit_name_field and unit_name_field != 'UNITNAME': |
347 | 352 | if unit_name_field in geology_gdf.columns: |
@@ -371,29 +376,36 @@ def calculate_thickness( |
371 | 376 | # Run thickness calculator |
372 | 377 | if calculator_type == "InterpolatedStructure": |
373 | 378 | calculator = InterpolatedStructure( |
374 | | - geology_gdf, |
375 | | - basal_contacts_gdf, |
376 | | - sampled_contacts_gdf, |
377 | | - structure_gdf, |
| 379 | + bounding_box=bounding_box, |
378 | 380 | dtm_data=dtm_gdal, |
379 | | - stratigraphic_column=stratigraphic_order, |
| 381 | + is_strike=orientation_type == 'Strike', |
| 382 | + max_line_length=max_line_length, |
380 | 383 | ) |
381 | 384 | else: # StructuralPoint |
382 | 385 | if max_line_length is None: |
383 | 386 | raise ValueError("max_line_length parameter is required for StructuralPoint calculator") |
384 | 387 | calculator = StructuralPoint( |
385 | | - geology_gdf, |
386 | | - basal_contacts_gdf, |
387 | | - sampled_contacts_gdf, |
388 | | - structure_gdf, |
389 | | - max_line_length=max_line_length, |
| 388 | + bounding_box=bounding_box, |
390 | 389 | dtm_data=dtm_gdal, |
391 | | - stratigraphic_column=stratigraphic_order, |
| 390 | + is_strike=orientation_type == 'Strike', |
| 391 | + max_line_length=max_line_length, |
392 | 392 | ) |
393 | | - |
394 | | - thickness = calculator.calculate_thickness() |
395 | | - |
| 393 | + if unit_name_field != 'UNITNAME' and unit_name_field in geology_gdf.columns: |
| 394 | + geology_gdf = geology_gdf.rename(columns={unit_name_field: 'UNITNAME'}) |
| 395 | + units = geology_gdf.copy() |
| 396 | + |
| 397 | + units_unique = units.drop_duplicates(subset=[unit_name_field]).reset_index(drop=True) |
| 398 | + units = pd.DataFrame({'name': units_unique[unit_name_field]}) |
| 399 | + basal_contacts_gdf['type'] = 'BASAL' # required by calculator |
| 400 | + thickness = calculator.compute( |
| 401 | + units, |
| 402 | + stratigraphic_order, |
| 403 | + basal_contacts_gdf, |
| 404 | + structure_gdf, |
| 405 | + geology_gdf, |
| 406 | + sampled_contacts_gdf, |
| 407 | + ) |
396 | 408 | if updater: |
397 | 409 | updater(f"Thickness calculation complete: {len(thickness)} records") |
398 | | - |
| 410 | + |
399 | 411 | return thickness |
0 commit comments