From c0bb5f27c48ab639017ab30ee65c00c2972f56a2 Mon Sep 17 00:00:00 2001 From: Lachlan Grose Date: Tue, 17 Jun 2025 13:09:35 +1000 Subject: [PATCH 01/20] template for basal contacts processing tool --- map2loop/processing/algorithms/__init__.py | 1 + .../algorithms/extract_basal_contacts.py | 76 +++++++++++++++++++ map2loop/processing/provider.py | 3 + 3 files changed, 80 insertions(+) create mode 100644 map2loop/processing/algorithms/__init__.py create mode 100644 map2loop/processing/algorithms/extract_basal_contacts.py diff --git a/map2loop/processing/algorithms/__init__.py b/map2loop/processing/algorithms/__init__.py new file mode 100644 index 0000000..618a57d --- /dev/null +++ b/map2loop/processing/algorithms/__init__.py @@ -0,0 +1 @@ +from .extract_basal_contacts import BasalContactsAlgorithm diff --git a/map2loop/processing/algorithms/extract_basal_contacts.py b/map2loop/processing/algorithms/extract_basal_contacts.py new file mode 100644 index 0000000..11d406c --- /dev/null +++ b/map2loop/processing/algorithms/extract_basal_contacts.py @@ -0,0 +1,76 @@ +""" +*************************************************************************** +* * +* This program is free software; you can redistribute it and/or modify * +* it under the terms of the GNU General Public License as published by * +* the Free Software Foundation; either version 2 of the License, or * +* (at your option) any later version. * +* * +*************************************************************************** +""" + +from typing import Any, Optional + +from qgis import processing +from qgis.core import ( + QgsFeatureSink, + QgsProcessing, + QgsProcessingAlgorithm, + QgsProcessingContext, + QgsProcessingException, + QgsProcessingFeedback, + QgsProcessingParameterFeatureSink, + QgsProcessingParameterFeatureSource, +) + + +class BasalContactsAlgorithm(QgsProcessingAlgorithm): + """Processing algorithm to create basal contacts.""" + + INPUT = "INPUT" + OUTPUT = "OUTPUT" + + def name(self) -> str: + """Return the algorithm name.""" + return "loop: basal_contacts" + + def displayName(self) -> str: + """Return the algorithm display name.""" + return "Loop3d: Basal Contacts" + + def group(self) -> str: + """Return the algorithm group name.""" + return "Loop3d" + + def groupId(self) -> str: + """Return the algorithm group ID.""" + return "loop3d" + + def initAlgorithm(self, config: Optional[dict[str, Any]] = None) -> None: + """Initialize the algorithm parameters.""" + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT, + "Geology Polygons", + [QgsProcessing.TypeVectorPolygon], + ) + ) + self.addParameter( + QgsProcessingParameterFeatureSink( + self.OUTPUT, + "Basal Contacts", + ) + ) + pass + + def processAlgorithm( + self, + parameters: dict[str, Any], + context: QgsProcessingContext, + feedback: QgsProcessingFeedback, + ) -> dict[str, Any]: + pass + + def createInstance(self) -> QgsProcessingAlgorithm: + """Create a new instance of the algorithm.""" + return self.__class__() # BasalContactsAlgorithm() diff --git a/map2loop/processing/provider.py b/map2loop/processing/provider.py index 6e0add8..ba879bb 100644 --- a/map2loop/processing/provider.py +++ b/map2loop/processing/provider.py @@ -16,6 +16,8 @@ __version__, ) +from .algorithms import BasalContactsAlgorithm + # ############################################################################ # ########## Classes ############### # ################################## @@ -26,6 +28,7 @@ class Map2LoopProvider(QgsProcessingProvider): def loadAlgorithms(self): """Loads all algorithms belonging to this provider.""" + self.addAlgorithm(BasalContactsAlgorithm()) pass def id(self) -> str: From f9878b0396c19f57160fc4c27b4086221e83c3b6 Mon Sep 17 00:00:00 2001 From: Lachlan Grose Date: Tue, 17 Jun 2025 13:11:57 +1000 Subject: [PATCH 02/20] Update linter.yml --- .github/workflows/linter.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/linter.yml b/.github/workflows/linter.yml index 57c2289..f87c4e5 100644 --- a/.github/workflows/linter.yml +++ b/.github/workflows/linter.yml @@ -7,7 +7,7 @@ on: pull_request: branches: - - master + - main paths: - '**.py' workflow_dispatch: @@ -55,5 +55,5 @@ jobs: token: ${{ secrets.GITHUB_TOKEN }} title: "style: auto format fixes" body: "This PR applies style fixes by black and ruff." - base: master + base: main branch: lint/style-fixes-${{ github.run_id }} From 773ef4844f12f3183da8caed672eb7af34dfbd35 Mon Sep 17 00:00:00 2001 From: rabii-chaarani Date: Mon, 23 Jun 2025 13:55:03 +0930 Subject: [PATCH 03/20] feature: Implement StratigraphySorterAlgorithm --- map2loop/processing/algorithms/sorter.py | 201 +++++++++++++++++++++++ 1 file changed, 201 insertions(+) create mode 100644 map2loop/processing/algorithms/sorter.py diff --git a/map2loop/processing/algorithms/sorter.py b/map2loop/processing/algorithms/sorter.py new file mode 100644 index 0000000..999fcad --- /dev/null +++ b/map2loop/processing/algorithms/sorter.py @@ -0,0 +1,201 @@ +from typing import Any, Optional + +from qgis import processing +from qgis.core import ( + QgsFeatureSink, + QgsFields, QgsField, QgsFeature, QgsGeometry, + QgsProcessing, + QgsProcessingAlgorithm, + QgsProcessingContext, + QgsProcessingException, + QgsProcessingFeedback, + QgsProcessingParameterEnum, + QgsProcessingParameterFeatureSink, + QgsProcessingParameterFeatureSource, + QgsVectorLayer, +) + +# ──────────────────────────────────────────────── +# map2loop sorters +# ──────────────────────────────────────────────── +from map2loop.map2loop.sorter import ( + SorterAlpha, + SorterAgeBased, + SorterMaximiseContacts, + SorterObservationProjections, + SorterUseNetworkX, + SorterUseHint, # kept for backwards compatibility +) + +# a lookup so we don’t need a giant if/else block +SORTER_LIST = { + "Age‐based": SorterAgeBased, + "NetworkX topological": SorterUseNetworkX, + "Hint (deprecated)": SorterUseHint, + "Adjacency α": SorterAlpha, + "Maximise contacts": SorterMaximiseContacts, + "Observation projections": SorterObservationProjections, +} + +class StratigraphySorterAlgorithm(QgsProcessingAlgorithm): + """ + Creates a one-column ‘stratigraphic column’ table ordered + by the selected map2loop sorter. + """ + + INPUT = "INPUT" + ALGO = "SORT_ALGO" + OUTPUT = "OUTPUT" + + # ---------------------------------------------------------- + # Metadata + # ---------------------------------------------------------- + def name(self) -> str: + return "loop_sorter" + + def displayName(self) -> str: + return "loop: Stratigraphic sorter" + + def group(self) -> str: + return "Loop3d" + + def groupId(self) -> str: + return "loop3d" + + # ---------------------------------------------------------- + # Parameters + # ---------------------------------------------------------- + def initAlgorithm(self, config: Optional[dict[str, Any]] = None) -> None: + + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT, + self.tr("Geology polygons"), + [QgsProcessing.TypeVectorPolygon], + ) + ) + + # enum so the user can pick the strategy from a dropdown + self.addParameter( + QgsProcessingParameterEnum( + self.ALGO, + self.tr("Sorting strategy"), + options=list(SORTER_LIST.keys()), + defaultValue=0, # Age-based is safest default + ) + ) #:contentReference[oaicite:0]{index=0} + + self.addParameter( + QgsProcessingParameterFeatureSink( + self.OUTPUT, + self.tr("Stratigraphic column"), + ) + ) + + # ---------------------------------------------------------- + # Core + # ---------------------------------------------------------- + def processAlgorithm( + self, + parameters: dict[str, Any], + context: QgsProcessingContext, + feedback: QgsProcessingFeedback, + ) -> dict[str, Any]: + + # 1 ► fetch user selections + in_layer: QgsVectorLayer = self.parameterAsVectorLayer(parameters, self.INPUT, context) + algo_index: int = self.parameterAsEnum(parameters, self.ALGO, context) + sorter_cls = list(SORTER_LIST.values())[algo_index] + + feedback.pushInfo(f"Using sorter: {sorter_cls.__name__}") + + # 2 ► convert QGIS layers / tables to pandas + # -------------------------------------------------- + # You must supply these three DataFrames: + # units_df — required (layerId, name, minAge, maxAge, group) + # relationships_df — required (Index1 / Unitname1, Index2 / Unitname2 …) + # contacts_df — required for all but Age‐based + # + # Typical workflow: + # • iterate over in_layer.getFeatures() + # • build dicts/lists + # • pd.DataFrame(…) + # + # NB: map2loop does *not* need geometries – only attribute values. + # -------------------------------------------------- + units_df, relationships_df, contacts_df, map_data = build_input_frames(in_layer, feedback) + + # 3 ► run the sorter + sorter = sorter_cls() # instantiation is always zero-argument + order = sorter.sort( + units_df, + relationships_df, + contacts_df, + map_data, + ) + + # 4 ► write an in-memory table with the result + sink_fields = QgsFields() + sink_fields.append(QgsField("strat_pos", int)) + sink_fields.append(QgsField("unit_name", str)) + + (sink, dest_id) = self.parameterAsSink( + parameters, + self.OUTPUT, + context, + sink_fields, + QgsWkbTypes.NoGeometry, + in_layer.sourceCrs(), + ) + + for pos, name in enumerate(order, start=1): + f = QgsFeature(sink_fields) + f.setAttributes([pos, name]) + sink.addFeature(f, QgsFeatureSink.FastInsert) + + return {self.OUTPUT: dest_id} + + # ---------------------------------------------------------- + def createInstance(self) -> QgsProcessingAlgorithm: + return StratigraphySorterAlgorithm() + + +# ------------------------------------------------------------------------- +# Helper stub – you must replace with *your* conversion logic +# ------------------------------------------------------------------------- +def build_input_frames(layer: QgsVectorLayer, feedback) -> tuple: + """ + Placeholder that turns the geology layer (and any other project + layers) into the four objects required by the sorter. + + Returns + ------- + (units_df, relationships_df, contacts_df, map_data) + """ + import pandas as pd + from map2loop.map2loop.mapdata import MapData # adjust import path if needed + + # Example: convert the geology layer to a very small units_df + units_records = [] + for f in layer.getFeatures(): + units_records.append( + dict( + layerId=f.id(), + name=f["UNITNAME"], # attribute names → your schema + minAge=f.attribute("MIN_AGE"), + maxAge=f.attribute("MAX_AGE"), + group=f["GROUP"], + ) + ) + units_df = pd.DataFrame.from_records(units_records) + + # relationships_df and contacts_df are domain-specific ─ fill them here + relationships_df = pd.DataFrame(columns=["Index1", "UNITNAME_1", "Index2", "UNITNAME_2"]) + contacts_df = pd.DataFrame(columns=["UNITNAME_1", "UNITNAME_2", "length"]) + + # map_data can be mocked if you only use Age-based sorter + map_data = MapData() # or MapData.from_project(…) / MapData.from_files(…) + + feedback.pushInfo(f"Units → {len(units_df)} records") + + return units_df, relationships_df, contacts_df, map_data From e7bb9f660416ef3a6ea9fa08f5a2df9a1fb91267 Mon Sep 17 00:00:00 2001 From: rabii-chaarani Date: Mon, 18 Aug 2025 12:35:07 +0930 Subject: [PATCH 04/20] feature: functions to convert layers to GeoDataFrame and DataFrame --- map2loop/main/vectorLayerWrapper.py | 78 +++++++++++++++++++++++++++++ 1 file changed, 78 insertions(+) create mode 100644 map2loop/main/vectorLayerWrapper.py diff --git a/map2loop/main/vectorLayerWrapper.py b/map2loop/main/vectorLayerWrapper.py new file mode 100644 index 0000000..7e146db --- /dev/null +++ b/map2loop/main/vectorLayerWrapper.py @@ -0,0 +1,78 @@ +import pandas as pd +import geopandas as gpd +from qgis.core import QgsRaster, QgsWkbTypes + + +def qgsLayerToGeoDataFrame(layer) -> gpd.GeoDataFrame: + if layer is None: + return None + features = layer.getFeatures() + fields = layer.fields() + data = {'geometry': []} + for f in fields: + data[f.name()] = [] + for feature in features: + geom = feature.geometry() + if geom.isEmpty(): + continue + data['geometry'].append(geom) + for f in fields: + data[f.name()].append(feature[f.name()]) + return gpd.GeoDataFrame(data, crs=layer.crs().authid()) + + +def qgsLayerToDataFrame(layer, dtm) -> pd.DataFrame: + """Convert a vector layer to a pandas DataFrame + samples the geometry using either points or the vertices of the lines + + :param layer: _description_ + :type layer: _type_ + :param dtm: Digital Terrain Model to evaluate Z values + :type dtm: _type_ or None + :return: the dataframe object + :rtype: pd.DataFrame + """ + if layer is None: + return None + fields = layer.fields() + data = {} + data['X'] = [] + data['Y'] = [] + data['Z'] = [] + + for field in fields: + data[field.name()] = [] + for feature in layer.getFeatures(): + geom = feature.geometry() + points = [] + if geom.isMultipart(): + if geom.type() == QgsWkbTypes.PointGeometry: + points = geom.asMultiPoint() + elif geom.type() == QgsWkbTypes.LineGeometry: + for line in geom.asMultiPolyline(): + points.extend(line) + # points = geom.asMultiPolyline()[0] + else: + if geom.type() == QgsWkbTypes.PointGeometry: + points = [geom.asPoint()] + elif geom.type() == QgsWkbTypes.LineGeometry: + points = geom.asPolyline() + + for p in points: + data['X'].append(p.x()) + data['Y'].append(p.y()) + if dtm is not None: + # Replace with your coordinates + + # Extract the value at the point + z_value = dtm.dataProvider().identify(p, QgsRaster.IdentifyFormatValue) + if z_value.isValid(): + z_value = z_value.results()[1] + else: + z_value = -9999 + data['Z'].append(z_value) + if dtm is None: + data['Z'].append(0) + for field in fields: + data[field.name()].append(feature[field.name()]) + return pd.DataFrame(data) From 95ba89b4544eea809501c1844060b1d78a47b32d Mon Sep 17 00:00:00 2001 From: rabii-chaarani Date: Mon, 18 Aug 2025 13:21:09 +0930 Subject: [PATCH 05/20] feature: add geodataframe to qgsLayer conversion --- map2loop/main/vectorLayerWrapper.py | 122 +++++++++++++++++++++++++++- 1 file changed, 121 insertions(+), 1 deletion(-) diff --git a/map2loop/main/vectorLayerWrapper.py b/map2loop/main/vectorLayerWrapper.py index 7e146db..81af364 100644 --- a/map2loop/main/vectorLayerWrapper.py +++ b/map2loop/main/vectorLayerWrapper.py @@ -1,6 +1,20 @@ +from qgis.core import ( + QgsVectorLayer, + QgsFields, + QgsField, + QgsFeature, + QgsGeometry, + QgsWkbTypes, + QgsCoordinateReferenceSystem, + QgsProject, + QgsRaster + ) +from qgis.PyQt.QtCore import QVariant, QDateTime + +from shapely.geometry import Point, MultiPoint, LineString, MultiLineString, Polygon, MultiPolygon import pandas as pd import geopandas as gpd -from qgis.core import QgsRaster, QgsWkbTypes + def qgsLayerToGeoDataFrame(layer) -> gpd.GeoDataFrame: @@ -76,3 +90,109 @@ def qgsLayerToDataFrame(layer, dtm) -> pd.DataFrame: for field in fields: data[field.name()].append(feature[field.name()]) return pd.DataFrame(data) + +def gdf_to_qgis_layer(gdf, layer_name="from_gdf"): + """ + Convert a GeoPandas GeoDataFrame to a QGIS memory layer (QgsVectorLayer). + Keeps attributes and CRS. Works for Point/LineString/Polygon and their Multi*. + """ + + if gdf is None or gdf.empty: + raise ValueError("GeoDataFrame is empty") + + # --- infer geometry type from first non-empty geometry + def infer_wkb(geoms): + for g in geoms: + if g is None: + continue + if hasattr(g, "is_empty") and g.is_empty: + continue + if isinstance(g, MultiPoint): return QgsWkbTypes.MultiPoint + if isinstance(g, Point): return QgsWkbTypes.Point + if isinstance(g, MultiLineString): return QgsWkbTypes.MultiLineString + if isinstance(g, LineString): return QgsWkbTypes.LineString + if isinstance(g, MultiPolygon): return QgsWkbTypes.MultiPolygon + if isinstance(g, Polygon): return QgsWkbTypes.Polygon + raise ValueError("Could not infer geometry type (all geometries empty?)") + + wkb_type = infer_wkb(gdf.geometry) + + # --- build CRS + crs_qgis = QgsCoordinateReferenceSystem() + if gdf.crs is not None: + try: + crs_qgis = QgsCoordinateReferenceSystem.fromWkt(gdf.crs.to_wkt()) + except Exception: + epsg = gdf.crs.to_epsg() + if epsg: + crs_qgis = QgsCoordinateReferenceSystem.fromEpsgId(int(epsg)) + + geom_str = QgsWkbTypes.displayString(wkb_type) # e.g. "LineString" + uri = f"{geom_str}?crs={crs_qgis.authid()}" if crs_qgis.isValid() else geom_str + layer = QgsVectorLayer(uri, layer_name, "memory") + prov = layer.dataProvider() + + # --- fields: map pandas dtypes → QGIS + import numpy as np + fields = QgsFields() + for col in gdf.columns: + if col == gdf.geometry.name: + continue + dtype = gdf[col].dtype + if pd.api.types.is_integer_dtype(dtype): + qtype = QVariant.Int + elif pd.api.types.is_float_dtype(dtype): + qtype = QVariant.Double + elif pd.api.types.is_bool_dtype(dtype): + qtype = QVariant.Bool + elif pd.api.types.is_datetime64_any_dtype(dtype): + qtype = QVariant.DateTime + else: + qtype = QVariant.String + fields.append(QgsField(str(col), qtype)) + prov.addAttributes(list(fields)) + layer.updateFields() + + # --- features + feats = [] + non_geom_cols = [c for c in gdf.columns if c != gdf.geometry.name] + + for _, row in gdf.iterrows(): + geom = row[gdf.geometry.name] + if geom is None or (hasattr(geom, "is_empty") and geom.is_empty): + continue + + f = QgsFeature(fields) + + # attributes in declared order with type cleanup + attrs = [] + for col in non_geom_cols: + val = row[col] + # numpy scalar → python scalar + if isinstance(val, (np.generic,)): + try: + val = val.item() + except Exception: + pass + # pandas Timestamp → QDateTime (if column is datetime) + if pd.api.types.is_datetime64_any_dtype(gdf[col].dtype): + if pd.isna(val): + val = None + else: + val = QDateTime(val.to_pydatetime()) + attrs.append(val) + f.setAttributes(attrs) + + # geometry (shapely → QGIS) + try: + f.setGeometry(QgsGeometry.fromWkb(geom.wkb)) + except Exception: + f.setGeometry(QgsGeometry.fromWkt(geom.wkt)) + + feats.append(f) + + if feats: + prov.addFeatures(feats) + layer.updateExtents() + + return layer # optionally: QgsProject.instance().addMapLayer(layer) From f744d83bb484218d11c3583a8154ef5c19e3f1ac Mon Sep 17 00:00:00 2001 From: rabii-chaarani Date: Mon, 18 Aug 2025 13:23:52 +0930 Subject: [PATCH 06/20] refactor: update GeoDataFrameToQgsLayer --- map2loop/main/vectorLayerWrapper.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/map2loop/main/vectorLayerWrapper.py b/map2loop/main/vectorLayerWrapper.py index 81af364..b106fb3 100644 --- a/map2loop/main/vectorLayerWrapper.py +++ b/map2loop/main/vectorLayerWrapper.py @@ -91,7 +91,7 @@ def qgsLayerToDataFrame(layer, dtm) -> pd.DataFrame: data[field.name()].append(feature[field.name()]) return pd.DataFrame(data) -def gdf_to_qgis_layer(gdf, layer_name="from_gdf"): +def GeoDataFrameToQgsLayer(gdf, layer_name="from_gdf"): """ Convert a GeoPandas GeoDataFrame to a QGIS memory layer (QgsVectorLayer). Keeps attributes and CRS. Works for Point/LineString/Polygon and their Multi*. From ea67c2e2819416cb81a642f9a6312ae8a0afeee7 Mon Sep 17 00:00:00 2001 From: rabii-chaarani Date: Mon, 18 Aug 2025 13:37:22 +0930 Subject: [PATCH 07/20] fix: GeoDataFrameToQgsLayer to FeatureSink --- map2loop/main/vectorLayerWrapper.py | 221 +++++++++++++++++++--------- 1 file changed, 150 insertions(+), 71 deletions(-) diff --git a/map2loop/main/vectorLayerWrapper.py b/map2loop/main/vectorLayerWrapper.py index b106fb3..90d0625 100644 --- a/map2loop/main/vectorLayerWrapper.py +++ b/map2loop/main/vectorLayerWrapper.py @@ -91,91 +91,170 @@ def qgsLayerToDataFrame(layer, dtm) -> pd.DataFrame: data[field.name()].append(feature[field.name()]) return pd.DataFrame(data) -def GeoDataFrameToQgsLayer(gdf, layer_name="from_gdf"): +def GeoDataFrameToQgsLayer(qgs_algorithm, geodataframe, parameters, context, output_key, feedback=None): """ - Convert a GeoPandas GeoDataFrame to a QGIS memory layer (QgsVectorLayer). - Keeps attributes and CRS. Works for Point/LineString/Polygon and their Multi*. + Write a GeoPandas GeoDataFrame directly to a QGIS Processing FeatureSink. + + Parameters + ---------- + alg : QgsProcessingAlgorithm (self) + gdf : geopandas.GeoDataFrame + parameters : dict (from processAlgorithm) + context : QgsProcessingContext + output_key : str (e.g. self.OUTPUT) + feedback : QgsProcessingFeedback | None + + Returns + ------- + str : dest_id to return from processAlgorithm, e.g. { output_key: dest_id } """ + import pandas as pd + import numpy as np + from shapely.geometry import ( + Point, MultiPoint, LineString, MultiLineString, Polygon, MultiPolygon + ) + + from qgis.core import ( + QgsFields, QgsField, QgsFeature, QgsGeometry, + QgsWkbTypes, QgsCoordinateReferenceSystem, QgsFeatureSink + ) + from qgis.PyQt.QtCore import QVariant, QDateTime + + if feedback is None: + class _Dummy: + def pushInfo(self, *a, **k): pass + def reportError(self, *a, **k): pass + def setProgress(self, *a, **k): pass + def isCanceled(self): return False + feedback = _Dummy() + + if geodataframe is None: + raise ValueError("GeoDataFrame is None") + if geodataframe.empty: + feedback.pushInfo("Input GeoDataFrame is empty; creating empty output layer.") + + # --- infer WKB type (family, Multi, Z) + def _infer_wkb(series): + base = None + any_multi = False + has_z = False + for geom in series: + if geom is None: continue + if getattr(geom, "is_empty", False): continue + # multi? + if isinstance(geom, (MultiPoint, MultiLineString, MultiPolygon)): + any_multi = True + g0 = next(iter(getattr(geom, "geoms", [])), None) + gt = getattr(g0, "geom_type", None) or None + else: + gt = getattr(geom, "geom_type", None) + + # base family + if gt in ("Point", "LineString", "Polygon"): + base = gt + # z? + try: + has_z = has_z or bool(getattr(geom, "has_z", False)) + except Exception: + pass + if base: + break + + if base is None: + # default safely to LineString if everything is empty; adjust if you prefer Point/Polygon + base = "LineString" + + fam = { + "Point": QgsWkbTypes.Point, + "LineString": QgsWkbTypes.LineString, + "Polygon": QgsWkbTypes.Polygon, + }[base] - if gdf is None or gdf.empty: - raise ValueError("GeoDataFrame is empty") - - # --- infer geometry type from first non-empty geometry - def infer_wkb(geoms): - for g in geoms: - if g is None: - continue - if hasattr(g, "is_empty") and g.is_empty: - continue - if isinstance(g, MultiPoint): return QgsWkbTypes.MultiPoint - if isinstance(g, Point): return QgsWkbTypes.Point - if isinstance(g, MultiLineString): return QgsWkbTypes.MultiLineString - if isinstance(g, LineString): return QgsWkbTypes.LineString - if isinstance(g, MultiPolygon): return QgsWkbTypes.MultiPolygon - if isinstance(g, Polygon): return QgsWkbTypes.Polygon - raise ValueError("Could not infer geometry type (all geometries empty?)") - - wkb_type = infer_wkb(gdf.geometry) - - # --- build CRS - crs_qgis = QgsCoordinateReferenceSystem() - if gdf.crs is not None: + if any_multi: + fam = QgsWkbTypes.multiType(fam) + if has_z: + fam = QgsWkbTypes.addZ(fam) + return fam + + wkb_type = _infer_wkb(geodataframe.geometry) + + # --- build CRS from gdf.crs + crs = QgsCoordinateReferenceSystem() + if geodataframe.crs is not None: try: - crs_qgis = QgsCoordinateReferenceSystem.fromWkt(gdf.crs.to_wkt()) + crs = QgsCoordinateReferenceSystem.fromWkt(geodataframe.crs.to_wkt()) except Exception: - epsg = gdf.crs.to_epsg() - if epsg: - crs_qgis = QgsCoordinateReferenceSystem.fromEpsgId(int(epsg)) + try: + epsg = geodataframe.crs.to_epsg() + if epsg: + crs = QgsCoordinateReferenceSystem.fromEpsgId(int(epsg)) + except Exception: + pass - geom_str = QgsWkbTypes.displayString(wkb_type) # e.g. "LineString" - uri = f"{geom_str}?crs={crs_qgis.authid()}" if crs_qgis.isValid() else geom_str - layer = QgsVectorLayer(uri, layer_name, "memory") - prov = layer.dataProvider() - - # --- fields: map pandas dtypes → QGIS - import numpy as np + # --- build QGIS fields from pandas dtypes fields = QgsFields() - for col in gdf.columns: - if col == gdf.geometry.name: - continue - dtype = gdf[col].dtype + non_geom_cols = [c for c in geodataframe.columns if c != geodataframe.geometry.name] + + def _qvariant_type(dtype) -> QVariant.Type: if pd.api.types.is_integer_dtype(dtype): - qtype = QVariant.Int - elif pd.api.types.is_float_dtype(dtype): - qtype = QVariant.Double - elif pd.api.types.is_bool_dtype(dtype): - qtype = QVariant.Bool - elif pd.api.types.is_datetime64_any_dtype(dtype): - qtype = QVariant.DateTime - else: - qtype = QVariant.String - fields.append(QgsField(str(col), qtype)) - prov.addAttributes(list(fields)) - layer.updateFields() - - # --- features - feats = [] - non_geom_cols = [c for c in gdf.columns if c != gdf.geometry.name] - - for _, row in gdf.iterrows(): - geom = row[gdf.geometry.name] - if geom is None or (hasattr(geom, "is_empty") and geom.is_empty): + return QVariant.Int + if pd.api.types.is_float_dtype(dtype): + return QVariant.Double + if pd.api.types.is_bool_dtype(dtype): + return QVariant.Bool + if pd.api.types.is_datetime64_any_dtype(dtype): + return QVariant.DateTime + return QVariant.String + + for col in non_geom_cols: + fields.append(QgsField(str(col), _qvariant_type(geodataframe[col].dtype))) + + # --- create sink + sink, dest_id = qgs_algorithm.parameterAsSink( + parameters, + output_key, + context, + fields, + wkb_type, + crs, + ) + if sink is None: + from qgis.core import QgsProcessingException + raise QgsProcessingException("Could not create output sink") + + # --- write features + total = len(geodataframe.index) + is_multi_sink = QgsWkbTypes.isMultiType(wkb_type) + + for i, (_, row) in enumerate(geodataframe.iterrows()): + if feedback.isCanceled(): + break + + geom = row[geodataframe.geometry.name] + if geom is None or getattr(geom, "is_empty", False): continue + # promote single → multi if needed + if is_multi_sink: + if isinstance(geom, Point): + geom = MultiPoint([geom]) + elif isinstance(geom, LineString): + geom = MultiLineString([geom]) + elif isinstance(geom, Polygon): + geom = MultiPolygon([geom]) + f = QgsFeature(fields) - # attributes in declared order with type cleanup + # attributes in declared order attrs = [] for col in non_geom_cols: val = row[col] - # numpy scalar → python scalar - if isinstance(val, (np.generic,)): + if isinstance(val, np.generic): try: val = val.item() except Exception: pass - # pandas Timestamp → QDateTime (if column is datetime) - if pd.api.types.is_datetime64_any_dtype(gdf[col].dtype): + if pd.api.types.is_datetime64_any_dtype(geodataframe[col].dtype): if pd.isna(val): val = None else: @@ -189,10 +268,10 @@ def infer_wkb(geoms): except Exception: f.setGeometry(QgsGeometry.fromWkt(geom.wkt)) - feats.append(f) + sink.addFeature(f, QgsFeatureSink.FastInsert) + + if total: + feedback.setProgress(int(100.0 * (i + 1) / total)) - if feats: - prov.addFeatures(feats) - layer.updateExtents() + return dest_id - return layer # optionally: QgsProject.instance().addMapLayer(layer) From f11d7131477b1ee28b6ae7effdec951541b3fef0 Mon Sep 17 00:00:00 2001 From: rabii-chaarani Date: Mon, 18 Aug 2025 13:37:36 +0930 Subject: [PATCH 08/20] refactor: clean up imports in vectorLayerWrapper.py --- map2loop/main/vectorLayerWrapper.py | 21 +++++---------------- 1 file changed, 5 insertions(+), 16 deletions(-) diff --git a/map2loop/main/vectorLayerWrapper.py b/map2loop/main/vectorLayerWrapper.py index 90d0625..e58b1b2 100644 --- a/map2loop/main/vectorLayerWrapper.py +++ b/map2loop/main/vectorLayerWrapper.py @@ -1,19 +1,19 @@ from qgis.core import ( - QgsVectorLayer, + QgsRaster, QgsFields, - QgsField, - QgsFeature, + QgsField, + QgsFeature, QgsGeometry, QgsWkbTypes, QgsCoordinateReferenceSystem, - QgsProject, - QgsRaster + QgsFeatureSink ) from qgis.PyQt.QtCore import QVariant, QDateTime from shapely.geometry import Point, MultiPoint, LineString, MultiLineString, Polygon, MultiPolygon import pandas as pd import geopandas as gpd +import numpy as np @@ -108,17 +108,6 @@ def GeoDataFrameToQgsLayer(qgs_algorithm, geodataframe, parameters, context, out ------- str : dest_id to return from processAlgorithm, e.g. { output_key: dest_id } """ - import pandas as pd - import numpy as np - from shapely.geometry import ( - Point, MultiPoint, LineString, MultiLineString, Polygon, MultiPolygon - ) - - from qgis.core import ( - QgsFields, QgsField, QgsFeature, QgsGeometry, - QgsWkbTypes, QgsCoordinateReferenceSystem, QgsFeatureSink - ) - from qgis.PyQt.QtCore import QVariant, QDateTime if feedback is None: class _Dummy: From d718b962120d91b310eef219acc9b25c8a789ce5 Mon Sep 17 00:00:00 2001 From: rabii-chaarani Date: Mon, 18 Aug 2025 13:38:02 +0930 Subject: [PATCH 09/20] feature: add basal contacts extraction algorithm --- .../algorithms/extract_basal_contacts.py | 61 ++++++++++++++++--- 1 file changed, 52 insertions(+), 9 deletions(-) diff --git a/map2loop/processing/algorithms/extract_basal_contacts.py b/map2loop/processing/algorithms/extract_basal_contacts.py index 11d406c..534958e 100644 --- a/map2loop/processing/algorithms/extract_basal_contacts.py +++ b/map2loop/processing/algorithms/extract_basal_contacts.py @@ -8,9 +8,10 @@ * * *************************************************************************** """ - +# Python imports from typing import Any, Optional +# QGIS imports from qgis import processing from qgis.core import ( QgsFeatureSink, @@ -22,17 +23,23 @@ QgsProcessingParameterFeatureSink, QgsProcessingParameterFeatureSource, ) +# Internal imports +from ...main.vectorLayerWrapper import qgsLayerToGeoDataFrame, GeoDataFrameToQgsLayer +from map2loop import ContactExtractor class BasalContactsAlgorithm(QgsProcessingAlgorithm): """Processing algorithm to create basal contacts.""" - - INPUT = "INPUT" - OUTPUT = "OUTPUT" + + + INPUT_GEOLOGY = 'GEOLOGY' + INPUT_FAULTS = 'FAULTS' + INPUT_STRATI_COLUMN = 'STRATIGRAPHIC_COLUMN' + OUTPUT = "BASAL_CONTACTS" def name(self) -> str: """Return the algorithm name.""" - return "loop: basal_contacts" + return "basal_contacts" def displayName(self) -> str: """Return the algorithm display name.""" @@ -48,20 +55,37 @@ def groupId(self) -> str: def initAlgorithm(self, config: Optional[dict[str, Any]] = None) -> None: """Initialize the algorithm parameters.""" + self.addParameter( QgsProcessingParameterFeatureSource( - self.INPUT, - "Geology Polygons", + self.INPUT_GEOLOGY, + "GEOLOGY", [QgsProcessing.TypeVectorPolygon], ) ) + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_FAULTS, + "FAULTS", + [QgsProcessing.TypeVectorLine], + optional=True, + ) + ) + + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_STRATI_COLUMN, + "STRATIGRAPHIC_COLUMN", + [QgsProcessing.TypeVectorLine], + ) + ) + self.addParameter( QgsProcessingParameterFeatureSink( self.OUTPUT, "Basal Contacts", ) ) - pass def processAlgorithm( self, @@ -69,7 +93,26 @@ def processAlgorithm( context: QgsProcessingContext, feedback: QgsProcessingFeedback, ) -> dict[str, Any]: - pass + + geology = self.parameterAsSource(parameters, self.INPUT_GEOLOGY, context) + faults = self.parameterAsSource(parameters, self.INPUT_FAULTS, context) + strati_column = self.parameterAsSource(parameters, self.INPUT_STRATI_COLUMN, context) + + geology = qgsLayerToGeoDataFrame(geology) + faults = qgsLayerToGeoDataFrame(faults) if faults else None + + feedback.pushInfo("Extracting Basal Contacts...") + contact_extractor = ContactExtractor(geology, faults, feedback) + contact_extractor.extract_basal_contacts(strati_column) + + basal_contacts = GeoDataFrameToQgsLayer( + self, + contact_extractor.basal_contacts, + parameters=parameters, + context=context, + feedback=feedback, + ) + return {self.OUTPUT: basal_contacts} def createInstance(self) -> QgsProcessingAlgorithm: """Create a new instance of the algorithm.""" From 0c7e7688c50f10ba74704c539070ac3984a96be6 Mon Sep 17 00:00:00 2001 From: rabii-chaarani Date: Mon, 18 Aug 2025 14:55:57 +0930 Subject: [PATCH 10/20] refactor: add imports in sorter.py --- map2loop/processing/algorithms/sorter.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/map2loop/processing/algorithms/sorter.py b/map2loop/processing/algorithms/sorter.py index 999fcad..a159855 100644 --- a/map2loop/processing/algorithms/sorter.py +++ b/map2loop/processing/algorithms/sorter.py @@ -3,7 +3,10 @@ from qgis import processing from qgis.core import ( QgsFeatureSink, - QgsFields, QgsField, QgsFeature, QgsGeometry, + QgsFields, + QgsField, + QgsFeature, + QgsGeometry, QgsProcessing, QgsProcessingAlgorithm, QgsProcessingContext, @@ -13,6 +16,7 @@ QgsProcessingParameterFeatureSink, QgsProcessingParameterFeatureSource, QgsVectorLayer, + QgsWkbTypes ) # ──────────────────────────────────────────────── From aeed68b750cbb46634a646f0cf67ae56d16d43d3 Mon Sep 17 00:00:00 2001 From: rabii-chaarani Date: Mon, 25 Aug 2025 14:18:44 +0930 Subject: [PATCH 11/20] feature: implement sampler algorithm --- map2loop/processing/algorithms/sampler.py | 149 ++++++++++++++++++++++ 1 file changed, 149 insertions(+) create mode 100644 map2loop/processing/algorithms/sampler.py diff --git a/map2loop/processing/algorithms/sampler.py b/map2loop/processing/algorithms/sampler.py new file mode 100644 index 0000000..5962a2a --- /dev/null +++ b/map2loop/processing/algorithms/sampler.py @@ -0,0 +1,149 @@ +""" +*************************************************************************** +* * +* This program is free software; you can redistribute it and/or modify * +* it under the terms of the GNU General Public License as published by * +* the Free Software Foundation; either version 2 of the License, or * +* (at your option) any later version. * +* * +*************************************************************************** +""" +# Python imports +from typing import Any, Optional + +# QGIS imports +from qgis import processing +from qgis.core import ( + QgsFeatureSink, + QgsProcessing, + QgsProcessingAlgorithm, + QgsProcessingContext, + QgsProcessingException, + QgsProcessingFeedback, + QgsProcessingParameterFeatureSink, + QgsProcessingParameterFeatureSource, + QgsProcessingParameterString, + QgsProcessingParameterNumber +) +# Internal imports +from ...main.vectorLayerWrapper import qgsLayerToGeoDataFrame, GeoDataFrameToQgsLayer +from map2loop.map2loop.sampler import SamplerDecimator, SamplerSpacing + + +class SamplerAlgorithm(QgsProcessingAlgorithm): + """Processing algorithm for sampling.""" + + INPUT_SAMPLER_TYPE = 'SAMPLER_TYPE' + INPUT_DTM = 'DTM' + INPUT_GEOLOGY = 'GEOLOGY' + INPUT_SPATIAL_DATA = 'SPATIAL_DATA' + INPUT_DECIMATION = 'DECIMATION' + INPUT_SPACING = 'SPACING' + + OUTPUT = "SAMPLED_CONTACTS" + + def name(self) -> str: + """Return the algorithm name.""" + return "sampler" + + def displayName(self) -> str: + """Return the algorithm display name.""" + return "Loop3d: Sampler" + + def group(self) -> str: + """Return the algorithm group name.""" + return "Loop3d" + + def groupId(self) -> str: + """Return the algorithm group ID.""" + return "loop3d" + + def initAlgorithm(self, config: Optional[dict[str, Any]] = None) -> None: + """Initialize the algorithm parameters.""" + + + self.addParameter( + QgsProcessingParameterString( + self.INPUT_SAMPLER_TYPE, + "SAMPLER_TYPE", + ) + ) + + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_DTM, + "DTM", + [QgsProcessing.TypeVectorRaster], + ) + ) + + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_GEOLOGY, + "GEOLOGY", + [QgsProcessing.TypeVectorPolygon], + optional=True, + ) + ) + + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_SPATIAL_DATA, + "SPATIAL_DATA", + [QgsProcessing.TypeVectorAnyGeometry], + optional=True, + ) + ) + + self.addParameter( + QgsProcessingParameterNumber( + self.INPUT_DECIMATION, + "DECIMATION", + optional=True, + ) + ) + + self.addParameter( + QgsProcessingParameterNumber( + self.INPUT_SPACING, + "SPACING", + optional=True, + ) + ) + + self.addParameter( + QgsProcessingParameterFeatureSink( + self.OUTPUT, + "Sampled Contacts", + ) + ) + + def processAlgorithm( + self, + parameters: dict[str, Any], + context: QgsProcessingContext, + feedback: QgsProcessingFeedback, + ) -> dict[str, Any]: + + dtm = self.parameterAsSource(parameters, self.INPUT_DTM, context) + geology = self.parameterAsSource(parameters, self.INPUT_GEOLOGY, context) + spatial_data = self.parameterAsSource(parameters, self.INPUT_SPATIAL_DATA, context) + decimation = self.parameterAsSource(parameters, self.INPUT_DECIMATION, context) + spacing = self.parameterAsSource(parameters, self.INPUT_SPACING, context) + sampler_type = self.parameterAsString(parameters, self.INPUT_SAMPLER_TYPE, context) + + # Convert geology layers to GeoDataFrames + geology = qgsLayerToGeoDataFrame(geology) + spatial_data = qgsLayerToGeoDataFrame(spatial_data) + + if sampler_type == "SamplerDecimator": + feedback.pushInfo("Sampling...") + sampler = SamplerDecimator(decimation=decimation, dtm_data=dtm, geology_data=geology, feedback=feedback) + samples = sampler.sample(spatial_data) + + samples = qgs + return {self.OUTPUT: basal_contacts} + + def createInstance(self) -> QgsProcessingAlgorithm: + """Create a new instance of the algorithm.""" + return self.__class__() # BasalContactsAlgorithm() \ No newline at end of file From 78f151123c69a48df9216dd21eb28c77d675638e Mon Sep 17 00:00:00 2001 From: rabii-chaarani Date: Mon, 25 Aug 2025 14:19:05 +0930 Subject: [PATCH 12/20] feature: add thickness calculator algorithms --- .../algorithms/thickness_calculator.py | 318 ++++++++++++++++++ 1 file changed, 318 insertions(+) create mode 100644 map2loop/processing/algorithms/thickness_calculator.py diff --git a/map2loop/processing/algorithms/thickness_calculator.py b/map2loop/processing/algorithms/thickness_calculator.py new file mode 100644 index 0000000..9b790ce --- /dev/null +++ b/map2loop/processing/algorithms/thickness_calculator.py @@ -0,0 +1,318 @@ +""" +*************************************************************************** +* * +* This program is free software; you can redistribute it and/or modify * +* it under the terms of the GNU General Public License as published by * +* the Free Software Foundation; either version 2 of the License, or * +* (at your option) any later version. * +* * +*************************************************************************** +""" +# Python imports +from typing import Any, Optional + +# QGIS imports +from qgis import processing +from qgis.core import ( + QgsFeatureSink, + QgsProcessing, + QgsProcessingAlgorithm, + QgsProcessingContext, + QgsProcessingException, + QgsProcessingFeedback, + QgsProcessingParameterFeatureSink, + QgsProcessingParameterFeatureSource, +) +# Internal imports +from ...main.vectorLayerWrapper import qgsLayerToGeoDataFrame, GeoDataFrameToQgsLayer +from map2loop.map2loop.thickness_calculator import InterpolatedStructure, StructuralPoint + + +class ThicknessCalculatorAlgorithm(QgsProcessingAlgorithm): + """Processing algorithm for thickness calculations.""" + + + INPUT_DTM = 'DTM' + INPUT_BOUNDING_BOX = 'BOUNDING_BOX' + INPUT_MAX_LINE_LENGTH = 'MAX_LINE_LENGTH' + INPUT_UNITS = 'UNITS' + INPUT_STRATI_COLUMN = 'STRATIGRAPHIC_COLUMN' + INPUT_BASAL_CONTACTS = 'BASAL_CONTACTS' + INPUT_STRUCTURE_DATA = 'STRUCTURE_DATA' + INPUT_GEOLOGY = 'GEOLOGY' + INPUT_SAMPLED_CONTACTS = 'SAMPLED_CONTACTS' + + OUTPUT = "THICKNESS" + + def name(self) -> str: + """Return the algorithm name.""" + return "thickness_calculator" + + def displayName(self) -> str: + """Return the algorithm display name.""" + return "Loop3d: Thickness Calculator" + + def group(self) -> str: + """Return the algorithm group name.""" + return "Loop3d" + + def groupId(self) -> str: + """Return the algorithm group ID.""" + return "loop3d" + + def initAlgorithm(self, config: Optional[dict[str, Any]] = None) -> None: + """Initialize the algorithm parameters.""" + + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_GEOLOGY, + "GEOLOGY", + [QgsProcessing.TypeVectorPolygon], + ) + ) + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_FAULTS, + "FAULTS", + [QgsProcessing.TypeVectorLine], + optional=True, + ) + ) + + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_STRATI_COLUMN, + "STRATIGRAPHIC_COLUMN", + [QgsProcessing.TypeVectorLine], + ) + ) + + self.addParameter( + QgsProcessingParameterFeatureSink( + self.OUTPUT, + "Basal Contacts", + ) + ) + + def processAlgorithm( + self, + parameters: dict[str, Any], + context: QgsProcessingContext, + feedback: QgsProcessingFeedback, + ) -> dict[str, Any]: + + geology = self.parameterAsSource(parameters, self.INPUT_GEOLOGY, context) + faults = self.parameterAsSource(parameters, self.INPUT_FAULTS, context) + strati_column = self.parameterAsSource(parameters, self.INPUT_STRATI_COLUMN, context) + + geology = qgsLayerToGeoDataFrame(geology) + faults = qgsLayerToGeoDataFrame(faults) if faults else None + + feedback.pushInfo("Extracting Basal Contacts...") + contact_extractor = ContactExtractor(geology, faults, feedback) + contact_extractor.extract_basal_contacts(strati_column) + + basal_contacts = GeoDataFrameToQgsLayer( + self, + contact_extractor.basal_contacts, + parameters=parameters, + context=context, + feedback=feedback, + ) + return {self.OUTPUT: basal_contacts} + + def createInstance(self) -> QgsProcessingAlgorithm: + """Create a new instance of the algorithm.""" + return self.__class__() # BasalContactsAlgorithm() + + + +class InterpolatedStructureAlgorithm(QgsProcessingAlgorithm): + """Processing algorithm for thickness calculations.""" + + + INPUT_UNITS = 'UNITS' + INPUT_STRATI_COLUMN = 'STRATIGRAPHIC_COLUMN' + INPUT_BASAL_CONTACTS = 'BASAL_CONTACTS' + INPUT_STRUCTURE_DATA = 'STRUCTURE_DATA' + INPUT_GEOLOGY = 'GEOLOGY' + INPUT_SAMPLED_CONTACTS = 'SAMPLED_CONTACTS' + + OUTPUT = "THICKNESS" + + def name(self) -> str: + """Return the algorithm name.""" + return "thickness_calculator" + + def displayName(self) -> str: + """Return the algorithm display name.""" + return "Loop3d: Thickness Calculator" + + def group(self) -> str: + """Return the algorithm group name.""" + return "Loop3d" + + def groupId(self) -> str: + """Return the algorithm group ID.""" + return "loop3d" + + def initAlgorithm(self, config: Optional[dict[str, Any]] = None) -> None: + """Initialize the algorithm parameters.""" + + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_GEOLOGY, + "GEOLOGY", + [QgsProcessing.TypeVectorPolygon], + ) + ) + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_FAULTS, + "FAULTS", + [QgsProcessing.TypeVectorLine], + optional=True, + ) + ) + + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_STRATI_COLUMN, + "STRATIGRAPHIC_COLUMN", + [QgsProcessing.TypeVectorLine], + ) + ) + + self.addParameter( + QgsProcessingParameterFeatureSink( + self.OUTPUT, + "Basal Contacts", + ) + ) + + def processAlgorithm( + self, + parameters: dict[str, Any], + context: QgsProcessingContext, + feedback: QgsProcessingFeedback, + ) -> dict[str, Any]: + + geology = self.parameterAsSource(parameters, self.INPUT_GEOLOGY, context) + faults = self.parameterAsSource(parameters, self.INPUT_FAULTS, context) + strati_column = self.parameterAsSource(parameters, self.INPUT_STRATI_COLUMN, context) + + geology = qgsLayerToGeoDataFrame(geology) + faults = qgsLayerToGeoDataFrame(faults) if faults else None + + feedback.pushInfo("Extracting Basal Contacts...") + contact_extractor = ContactExtractor(geology, faults, feedback) + contact_extractor.extract_basal_contacts(strati_column) + + basal_contacts = GeoDataFrameToQgsLayer( + self, + contact_extractor.basal_contacts, + parameters=parameters, + context=context, + feedback=feedback, + ) + return {self.OUTPUT: basal_contacts} + + def createInstance(self) -> QgsProcessingAlgorithm: + """Create a new instance of the algorithm.""" + return self.__class__() # BasalContactsAlgorithm() + + + +class StructuralPointAlgorithm(QgsProcessingAlgorithm): + """Processing algorithm for thickness calculations.""" + + + INPUT_UNITS = 'UNITS' + INPUT_STRATI_COLUMN = 'STRATIGRAPHIC_COLUMN' + INPUT_BASAL_CONTACTS = 'BASAL_CONTACTS' + INPUT_STRUCTURE_DATA = 'STRUCTURE_DATA' + INPUT_GEOLOGY = 'GEOLOGY' + INPUT_SAMPLED_CONTACTS = 'SAMPLED_CONTACTS' + + OUTPUT = "THICKNESS" + + def name(self) -> str: + """Return the algorithm name.""" + return "thickness_calculator" + + def displayName(self) -> str: + """Return the algorithm display name.""" + return "Loop3d: Thickness Calculator" + + def group(self) -> str: + """Return the algorithm group name.""" + return "Loop3d" + + def groupId(self) -> str: + """Return the algorithm group ID.""" + return "loop3d" + + def initAlgorithm(self, config: Optional[dict[str, Any]] = None) -> None: + """Initialize the algorithm parameters.""" + + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_GEOLOGY, + "GEOLOGY", + [QgsProcessing.TypeVectorPolygon], + ) + ) + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_FAULTS, + "FAULTS", + [QgsProcessing.TypeVectorLine], + optional=True, + ) + ) + + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_STRATI_COLUMN, + "STRATIGRAPHIC_COLUMN", + [QgsProcessing.TypeVectorLine], + ) + ) + + self.addParameter( + QgsProcessingParameterFeatureSink( + self.OUTPUT, + "Basal Contacts", + ) + ) + + def processAlgorithm( + self, + parameters: dict[str, Any], + context: QgsProcessingContext, + feedback: QgsProcessingFeedback, + ) -> dict[str, Any]: + + geology = self.parameterAsSource(parameters, self.INPUT_GEOLOGY, context) + faults = self.parameterAsSource(parameters, self.INPUT_FAULTS, context) + strati_column = self.parameterAsSource(parameters, self.INPUT_STRATI_COLUMN, context) + + geology = qgsLayerToGeoDataFrame(geology) + faults = qgsLayerToGeoDataFrame(faults) if faults else None + + feedback.pushInfo("Extracting Basal Contacts...") + contact_extractor = ContactExtractor(geology, faults, feedback) + contact_extractor.extract_basal_contacts(strati_column) + + basal_contacts = GeoDataFrameToQgsLayer( + self, + contact_extractor.basal_contacts, + parameters=parameters, + context=context, + feedback=feedback, + ) + return {self.OUTPUT: basal_contacts} + + def createInstance(self) -> QgsProcessingAlgorithm: + """Create a new instance of the algorithm.""" + return self.__class__() # BasalContactsAlgorithm() \ No newline at end of file From 40694065ebf7a0d40a6312f2e2f45aa1742c67a1 Mon Sep 17 00:00:00 2001 From: rabii-chaarani Date: Mon, 25 Aug 2025 14:25:55 +0930 Subject: [PATCH 13/20] fix: correct return value --- map2loop/processing/algorithms/sampler.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/map2loop/processing/algorithms/sampler.py b/map2loop/processing/algorithms/sampler.py index 5962a2a..be6b134 100644 --- a/map2loop/processing/algorithms/sampler.py +++ b/map2loop/processing/algorithms/sampler.py @@ -142,8 +142,8 @@ def processAlgorithm( samples = sampler.sample(spatial_data) samples = qgs - return {self.OUTPUT: basal_contacts} + return {self.OUTPUT: samples} def createInstance(self) -> QgsProcessingAlgorithm: """Create a new instance of the algorithm.""" - return self.__class__() # BasalContactsAlgorithm() \ No newline at end of file + return self.__class__() # SamplerAlgorithm() \ No newline at end of file From f67be330d9a63cb7543ab7e9b4e58e2e252e3cfa Mon Sep 17 00:00:00 2001 From: rabii-chaarani Date: Mon, 25 Aug 2025 14:27:31 +0930 Subject: [PATCH 14/20] feature: add support SamplerSpacing --- map2loop/processing/algorithms/sampler.py | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/map2loop/processing/algorithms/sampler.py b/map2loop/processing/algorithms/sampler.py index be6b134..b9a182b 100644 --- a/map2loop/processing/algorithms/sampler.py +++ b/map2loop/processing/algorithms/sampler.py @@ -137,11 +137,16 @@ def processAlgorithm( spatial_data = qgsLayerToGeoDataFrame(spatial_data) if sampler_type == "SamplerDecimator": - feedback.pushInfo("Sampling...") - sampler = SamplerDecimator(decimation=decimation, dtm_data=dtm, geology_data=geology, feedback=feedback) - samples = sampler.sample(spatial_data) - - samples = qgs + feedback.pushInfo("Sampling...") + sampler = SamplerDecimator(decimation=decimation, dtm_data=dtm, geology_data=geology, feedback=feedback) + samples = sampler.sample(spatial_data) + if sampler_type == "SamplerSpacing": + feedback.pushInfo("Sampling...") + sampler = SamplerSpacing(spacing=spacing, dtm_data=dtm, geology_data=geology, feedback=feedback) + samples = sampler.sample(spatial_data) + + #TODO: convert sample to qgis layer + # samples = qgs return {self.OUTPUT: samples} def createInstance(self) -> QgsProcessingAlgorithm: From 3c465534b83e1d61789e479d49850484e66b4ca6 Mon Sep 17 00:00:00 2001 From: rabii-chaarani Date: Mon, 25 Aug 2025 14:28:40 +0930 Subject: [PATCH 15/20] fix: update sampler type strings for consistency --- map2loop/processing/algorithms/sampler.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/map2loop/processing/algorithms/sampler.py b/map2loop/processing/algorithms/sampler.py index b9a182b..1011cf6 100644 --- a/map2loop/processing/algorithms/sampler.py +++ b/map2loop/processing/algorithms/sampler.py @@ -136,11 +136,12 @@ def processAlgorithm( geology = qgsLayerToGeoDataFrame(geology) spatial_data = qgsLayerToGeoDataFrame(spatial_data) - if sampler_type == "SamplerDecimator": + if sampler_type == "decimator": feedback.pushInfo("Sampling...") sampler = SamplerDecimator(decimation=decimation, dtm_data=dtm, geology_data=geology, feedback=feedback) samples = sampler.sample(spatial_data) - if sampler_type == "SamplerSpacing": + + if sampler_type == "spacing": feedback.pushInfo("Sampling...") sampler = SamplerSpacing(spacing=spacing, dtm_data=dtm, geology_data=geology, feedback=feedback) samples = sampler.sample(spatial_data) From 7e76267655033124d0434498d385306f4145ccb4 Mon Sep 17 00:00:00 2001 From: rabii-chaarani Date: Tue, 26 Aug 2025 11:47:41 +0930 Subject: [PATCH 16/20] fix: convert dataframe to qgis layer --- map2loop/processing/algorithms/sampler.py | 48 ++++++++++++++++++++--- 1 file changed, 42 insertions(+), 6 deletions(-) diff --git a/map2loop/processing/algorithms/sampler.py b/map2loop/processing/algorithms/sampler.py index 1011cf6..3feb3e6 100644 --- a/map2loop/processing/algorithms/sampler.py +++ b/map2loop/processing/algorithms/sampler.py @@ -10,9 +10,9 @@ """ # Python imports from typing import Any, Optional +from qgis.PyQt.QtCore import QMetaType # QGIS imports -from qgis import processing from qgis.core import ( QgsFeatureSink, QgsProcessing, @@ -23,10 +23,15 @@ QgsProcessingParameterFeatureSink, QgsProcessingParameterFeatureSource, QgsProcessingParameterString, - QgsProcessingParameterNumber + QgsProcessingParameterNumber, + QgsField, + QgsFeature, + QgsGeometry, + QgsPointXY, + QgsVectorLayer ) # Internal imports -from ...main.vectorLayerWrapper import qgsLayerToGeoDataFrame, GeoDataFrameToQgsLayer +from ...main.vectorLayerWrapper import qgsLayerToGeoDataFrame from map2loop.map2loop.sampler import SamplerDecimator, SamplerSpacing @@ -145,10 +150,41 @@ def processAlgorithm( feedback.pushInfo("Sampling...") sampler = SamplerSpacing(spacing=spacing, dtm_data=dtm, geology_data=geology, feedback=feedback) samples = sampler.sample(spatial_data) + - #TODO: convert sample to qgis layer - # samples = qgs - return {self.OUTPUT: samples} + # create layer + vector_layer = QgsVectorLayer("Point", "sampled_points", "memory") + provider = vector_layer.dataProvider() + + # add fields + provider.addAttributes([QgsField("ID", QMetaType.Type.QString), + QgsField("X", QMetaType.Type.Float), + QgsField("Y", QMetaType.Type.Float), + QgsField("Z", QMetaType.Type.Float), + QgsField("featureId", QMetaType.Type.QString) + ]) + vector_layer.updateFields() # tell the vector layer to fetch changes from the provider + + # add a feature + for i in range(len(samples)): + feature = QgsFeature() + feature.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(samples.X[i], samples.Y[i], samples.Z[i]))) + feature.setAttributes([samples.ID[i], samples.X[i], samples.Y[i], samples.Z[i], samples.featureId[i]]) + provider.addFeatures([feature]) + + # update layer's extent when new features have been added + # because change of extent in provider is not propagated to the layer + vector_layer.updateExtents() + # --- create sink + sink, dest_id = self.parameterAsSink( + parameters, + self.OUTPUT, + context, + vector_layer.fields(), + QgsGeometry.Type.Point, + spatial_data.crs, + ) + return {self.OUTPUT: dest_id} def createInstance(self) -> QgsProcessingAlgorithm: """Create a new instance of the algorithm.""" From c42a4ddfc86a9adb1191f9160de25d884193c502 Mon Sep 17 00:00:00 2001 From: rabii-chaarani Date: Wed, 27 Aug 2025 10:36:57 +0930 Subject: [PATCH 17/20] feature: add dataframe to point sink conversion --- map2loop/main/vectorLayerWrapper.py | 189 +++++++++++++++++++++++++++- 1 file changed, 187 insertions(+), 2 deletions(-) diff --git a/map2loop/main/vectorLayerWrapper.py b/map2loop/main/vectorLayerWrapper.py index e58b1b2..7069c88 100644 --- a/map2loop/main/vectorLayerWrapper.py +++ b/map2loop/main/vectorLayerWrapper.py @@ -1,3 +1,5 @@ +# PyQGIS / PyQt imports + from qgis.core import ( QgsRaster, QgsFields, @@ -6,8 +8,10 @@ QgsGeometry, QgsWkbTypes, QgsCoordinateReferenceSystem, - QgsFeatureSink + QgsFeatureSink, + QgsProcessingException ) + from qgis.PyQt.QtCore import QVariant, QDateTime from shapely.geometry import Point, MultiPoint, LineString, MultiLineString, Polygon, MultiPolygon @@ -208,7 +212,6 @@ def _qvariant_type(dtype) -> QVariant.Type: crs, ) if sink is None: - from qgis.core import QgsProcessingException raise QgsProcessingException("Could not create output sink") # --- write features @@ -264,3 +267,185 @@ def _qvariant_type(dtype) -> QVariant.Type: return dest_id + +# ---------- helpers ---------- + +def _qvariant_type_from_dtype(dtype) -> QVariant.Type: + """Map a pandas dtype to a QVariant type.""" + import numpy as np + if np.issubdtype(dtype, np.integer): + # prefer 64-bit when detected + try: + return QVariant.LongLong + except AttributeError: + return QVariant.Int + if np.issubdtype(dtype, np.floating): + return QVariant.Double + if np.issubdtype(dtype, np.bool_): + return QVariant.Bool + # datetimes + try: + import pandas as pd + if pd.api.types.is_datetime64_any_dtype(dtype): + return QVariant.DateTime + if pd.api.types.is_datetime64_ns_dtype(dtype): + return QVariant.DateTime + if pd.api.types.is_datetime64_dtype(dtype): + return QVariant.DateTime + if pd.api.types.is_timedelta64_dtype(dtype): + # store as string "HH:MM:SS" fallback + return QVariant.String + except Exception: + pass + # default to string + return QVariant.String + + +def _fields_from_dataframe(df, drop_cols=None) -> QgsFields: + """Build QgsFields from DataFrame dtypes.""" + drop_cols = set(drop_cols or []) + fields = QgsFields() + for name, dtype in df.dtypes.items(): + if name in drop_cols: + continue + vtype = _qvariant_type_from_dtype(dtype) + fields.append(QgsField(name, vtype)) + return fields + + +# ---------- main function you'll call inside processAlgorithm ---------- + +def dataframe_to_point_sink( + df, + x_col: str, + y_col: str, + *, + crs: QgsCoordinateReferenceSystem, + algorithm, # `self` inside a QgsProcessingAlgorithm + parameters: dict, + context, + feedback, + sink_param_name: str = "OUTPUT", + z_col: str = None, + m_col: str = None, + include_coords_in_attrs: bool = False, +): + """ + Write a pandas DataFrame to a point feature sink (QgsProcessingParameterFeatureSink). + + Params + ------ + df : pandas.DataFrame Data with coordinate columns. + x_col, y_col : str Column names for X/Easting/Longitude and Y/Northing/Latitude. + crs : QgsCoordinateReferenceSystem CRS of the coordinates (e.g., QgsCoordinateReferenceSystem('EPSG:4326')). + algorithm : QgsProcessingAlgorithm Use `self` from inside processAlgorithm. + parameters, context, feedback Standard Processing plumbing. + sink_param_name : str Name of your sink output parameter (default "OUTPUT"). + z_col, m_col : str | None Optional Z and M columns for 3D/M points. + include_coords_in_attrs : bool If False, x/y/z/m are not written as attributes. + + Returns + ------- + (sink, sink_id) The created sink and its ID. Also returns feature count via feedback. + """ + import pandas as pd + if not isinstance(df, pd.DataFrame): + raise TypeError("df must be a pandas.DataFrame") + + # Make a working copy; optionally drop coordinate columns from attributes + attr_df = df.copy() + drop_cols = [] + for col in [x_col, y_col, z_col, m_col]: + if col and not include_coords_in_attrs: + drop_cols.append(col) + + fields = _fields_from_dataframe(attr_df, drop_cols=drop_cols) + + # Geometry type (2D/3D/M) + has_z = z_col is not None and z_col in df.columns + has_m = m_col is not None and m_col in df.columns + if has_z and has_m: + wkb = QgsWkbTypes.PointZM + elif has_z: + wkb = QgsWkbTypes.PointZ + elif has_m: + wkb = QgsWkbTypes.PointM + else: + wkb = QgsWkbTypes.Point + + # Create the sink + sink, sink_id = algorithm.parameterAsSink( + parameters, + sink_param_name, + context, + fields, + wkb, + crs + ) + if sink is None: + raise QgsProcessingException("Could not create feature sink. Check output parameter and inputs.") + + total = len(df) + feedback.pushInfo(f"Writing {total} features…") + + # Precompute attribute column order + attr_columns = [f.name() for f in fields] + + # Iterate rows and write features + for i, (idx, row) in enumerate(df.iterrows(), start=1): + if feedback.isCanceled(): + break + + # Build point geometry + x = row[x_col] + y = row[y_col] + + # skip rows with missing coords + if pd.isna(x) or pd.isna(y): + continue + + if has_z and not pd.isna(row[z_col]) and has_m and not pd.isna(row[m_col]): + pt = QgsPoint(float(x), float(y), float(row[z_col]), float(row[m_col])) + elif has_z and not pd.isna(row[z_col]): + pt = QgsPoint(float(x), float(y), float(row[z_col])) + elif has_m and not pd.isna(row[m_col]): + # PointM constructor: setZValue not needed; M is the 4th ordinate + pt = QgsPoint(float(x), float(y)) + pt.setM(float(row[m_col])) + else: + pt = QgsPointXY(float(x), float(y)) + + feat = QgsFeature(fields) + feat.setGeometry(QgsGeometry.fromPoint(pt) if isinstance(pt, QgsPoint) else QgsGeometry.fromPointXY(pt)) + + # Attributes in the same order as fields + attrs = [] + for col in attr_columns: + val = row[col] if col in row else None + # Pandas NaN -> None + if pd.isna(val): + val = None + # Convert numpy types to Python scalars to avoid QVariant issues + try: + import numpy as np + if isinstance(val, (np.generic,)): + val = val.item() + except Exception: + pass + # Convert pandas Timestamp to Python datetime + if hasattr(val, "to_pydatetime"): + try: + val = val.to_pydatetime() + except Exception: + val = str(val) + attrs.append(val) + feat.setAttributes(attrs) + + sink.addFeature(feat, QgsFeature.FastInsert) + + if i % 1000 == 0: + feedback.setProgress(int(100.0 * i / max(total, 1))) + + feedback.pushInfo("Done.") + feedback.setProgress(100) + return sink, sink_id From 54a0bc98eea26f2b22632e4e8f8c03c66798c360 Mon Sep 17 00:00:00 2001 From: rabii-chaarani Date: Wed, 27 Aug 2025 10:37:27 +0930 Subject: [PATCH 18/20] fix: add input parameters --- .../algorithms/thickness_calculator.py | 30 +++++++++++++++++-- 1 file changed, 28 insertions(+), 2 deletions(-) diff --git a/map2loop/processing/algorithms/thickness_calculator.py b/map2loop/processing/algorithms/thickness_calculator.py index 9b790ce..24b684a 100644 --- a/map2loop/processing/algorithms/thickness_calculator.py +++ b/map2loop/processing/algorithms/thickness_calculator.py @@ -31,7 +31,7 @@ class ThicknessCalculatorAlgorithm(QgsProcessingAlgorithm): """Processing algorithm for thickness calculations.""" - + INPUT_THICKNESS_CALCULATOR_TYPE = 'THICKNESS_CALCULATOR_TYPE' INPUT_DTM = 'DTM' INPUT_BOUNDING_BOX = 'BOUNDING_BOX' INPUT_MAX_LINE_LENGTH = 'MAX_LINE_LENGTH' @@ -62,7 +62,33 @@ def groupId(self) -> str: def initAlgorithm(self, config: Optional[dict[str, Any]] = None) -> None: """Initialize the algorithm parameters.""" + + + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_THICKNESS_CALCULATOR_TYPE, + "Thickness Calculator Type", + [QgsProcessing.TypeVectorPoint], + ) + ) + + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_BASAL_CONTACTS, + "Basal Contacts", + [QgsProcessing.TypeVectorPoint], + ) + ) + + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_DTM, + "DTM", + [QgsProcessing.TypeVectorRaster], + ) + ) + self.addParameter( QgsProcessingParameterFeatureSource( self.INPUT_GEOLOGY, @@ -90,7 +116,7 @@ def initAlgorithm(self, config: Optional[dict[str, Any]] = None) -> None: self.addParameter( QgsProcessingParameterFeatureSink( self.OUTPUT, - "Basal Contacts", + "Thickness", ) ) From df54c1b16304bb386dc2ddb988fa4951a8e58f7f Mon Sep 17 00:00:00 2001 From: rabii-chaarani Date: Wed, 27 Aug 2025 11:58:30 +0930 Subject: [PATCH 19/20] refactor: rename to dataframeToQgsLayer --- map2loop/main/vectorLayerWrapper.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/map2loop/main/vectorLayerWrapper.py b/map2loop/main/vectorLayerWrapper.py index 7069c88..d68ed3d 100644 --- a/map2loop/main/vectorLayerWrapper.py +++ b/map2loop/main/vectorLayerWrapper.py @@ -315,7 +315,7 @@ def _fields_from_dataframe(df, drop_cols=None) -> QgsFields: # ---------- main function you'll call inside processAlgorithm ---------- -def dataframe_to_point_sink( +def dataframeToQgsLayer( df, x_col: str, y_col: str, From 703a89c6baf2561d0f4eee6406c43905b8ba4bed Mon Sep 17 00:00:00 2001 From: rabii-chaarani Date: Wed, 27 Aug 2025 11:58:42 +0930 Subject: [PATCH 20/20] refactor: update thickness calculator parameters and processing logic --- .../algorithms/thickness_calculator.py | 308 ++++++------------ 1 file changed, 93 insertions(+), 215 deletions(-) diff --git a/map2loop/processing/algorithms/thickness_calculator.py b/map2loop/processing/algorithms/thickness_calculator.py index 24b684a..f56b09c 100644 --- a/map2loop/processing/algorithms/thickness_calculator.py +++ b/map2loop/processing/algorithms/thickness_calculator.py @@ -22,9 +22,13 @@ QgsProcessingFeedback, QgsProcessingParameterFeatureSink, QgsProcessingParameterFeatureSource, + QgsProcessingParameterEnum, + QgsProcessingParameterNumber, + QgsProcessingParameterField, + QgsProcessingParameterMatrix ) # Internal imports -from ...main.vectorLayerWrapper import qgsLayerToGeoDataFrame, GeoDataFrameToQgsLayer +from ...main.vectorLayerWrapper import qgsLayerToGeoDataFrame, GeoDataFrameToQgsLayer, qgsLayerToDataFrame, dataframeToQgsLayer from map2loop.map2loop.thickness_calculator import InterpolatedStructure, StructuralPoint @@ -62,25 +66,15 @@ def groupId(self) -> str: def initAlgorithm(self, config: Optional[dict[str, Any]] = None) -> None: """Initialize the algorithm parameters.""" - - self.addParameter( - QgsProcessingParameterFeatureSource( + QgsProcessingParameterEnum( self.INPUT_THICKNESS_CALCULATOR_TYPE, "Thickness Calculator Type", - [QgsProcessing.TypeVectorPoint], + options=['InterpolatedStructure','StructuralPoint'], + allowMultiple=False, ) ) - - self.addParameter( - QgsProcessingParameterFeatureSource( - self.INPUT_BASAL_CONTACTS, - "Basal Contacts", - [QgsProcessing.TypeVectorPoint], - ) - ) - self.addParameter( QgsProcessingParameterFeatureSource( self.INPUT_DTM, @@ -88,103 +82,36 @@ def initAlgorithm(self, config: Optional[dict[str, Any]] = None) -> None: [QgsProcessing.TypeVectorRaster], ) ) - self.addParameter( - QgsProcessingParameterFeatureSource( - self.INPUT_GEOLOGY, - "GEOLOGY", - [QgsProcessing.TypeVectorPolygon], + QgsProcessingParameterEnum( + self.INPUT_BOUNDING_BOX, + "Bounding Box", + options=['minx','miny','maxx','maxy'], + allowMultiple=True, ) ) + self.addParameter( + QgsProcessingParameterNumber( + self.INPUT_MAX_LINE_LENGTH, + "Max Line Length", + minValue=0, + defaultValue=1000 + ) + ) self.addParameter( QgsProcessingParameterFeatureSource( - self.INPUT_FAULTS, - "FAULTS", + self.INPUT_UNITS, + "Units", [QgsProcessing.TypeVectorLine], - optional=True, ) ) - self.addParameter( QgsProcessingParameterFeatureSource( - self.INPUT_STRATI_COLUMN, - "STRATIGRAPHIC_COLUMN", + self.INPUT_BASAL_CONTACTS, + "Basal Contacts", [QgsProcessing.TypeVectorLine], ) ) - - self.addParameter( - QgsProcessingParameterFeatureSink( - self.OUTPUT, - "Thickness", - ) - ) - - def processAlgorithm( - self, - parameters: dict[str, Any], - context: QgsProcessingContext, - feedback: QgsProcessingFeedback, - ) -> dict[str, Any]: - - geology = self.parameterAsSource(parameters, self.INPUT_GEOLOGY, context) - faults = self.parameterAsSource(parameters, self.INPUT_FAULTS, context) - strati_column = self.parameterAsSource(parameters, self.INPUT_STRATI_COLUMN, context) - - geology = qgsLayerToGeoDataFrame(geology) - faults = qgsLayerToGeoDataFrame(faults) if faults else None - - feedback.pushInfo("Extracting Basal Contacts...") - contact_extractor = ContactExtractor(geology, faults, feedback) - contact_extractor.extract_basal_contacts(strati_column) - - basal_contacts = GeoDataFrameToQgsLayer( - self, - contact_extractor.basal_contacts, - parameters=parameters, - context=context, - feedback=feedback, - ) - return {self.OUTPUT: basal_contacts} - - def createInstance(self) -> QgsProcessingAlgorithm: - """Create a new instance of the algorithm.""" - return self.__class__() # BasalContactsAlgorithm() - - - -class InterpolatedStructureAlgorithm(QgsProcessingAlgorithm): - """Processing algorithm for thickness calculations.""" - - - INPUT_UNITS = 'UNITS' - INPUT_STRATI_COLUMN = 'STRATIGRAPHIC_COLUMN' - INPUT_BASAL_CONTACTS = 'BASAL_CONTACTS' - INPUT_STRUCTURE_DATA = 'STRUCTURE_DATA' - INPUT_GEOLOGY = 'GEOLOGY' - INPUT_SAMPLED_CONTACTS = 'SAMPLED_CONTACTS' - - OUTPUT = "THICKNESS" - - def name(self) -> str: - """Return the algorithm name.""" - return "thickness_calculator" - - def displayName(self) -> str: - """Return the algorithm display name.""" - return "Loop3d: Thickness Calculator" - - def group(self) -> str: - """Return the algorithm group name.""" - return "Loop3d" - - def groupId(self) -> str: - """Return the algorithm group ID.""" - return "loop3d" - - def initAlgorithm(self, config: Optional[dict[str, Any]] = None) -> None: - """Initialize the algorithm parameters.""" - self.addParameter( QgsProcessingParameterFeatureSource( self.INPUT_GEOLOGY, @@ -193,122 +120,32 @@ def initAlgorithm(self, config: Optional[dict[str, Any]] = None) -> None: ) ) self.addParameter( - QgsProcessingParameterFeatureSource( - self.INPUT_FAULTS, - "FAULTS", - [QgsProcessing.TypeVectorLine], - optional=True, - ) - ) - - self.addParameter( - QgsProcessingParameterFeatureSource( - self.INPUT_STRATI_COLUMN, - "STRATIGRAPHIC_COLUMN", - [QgsProcessing.TypeVectorLine], - ) - ) - - self.addParameter( - QgsProcessingParameterFeatureSink( - self.OUTPUT, - "Basal Contacts", - ) - ) - - def processAlgorithm( - self, - parameters: dict[str, Any], - context: QgsProcessingContext, - feedback: QgsProcessingFeedback, - ) -> dict[str, Any]: - - geology = self.parameterAsSource(parameters, self.INPUT_GEOLOGY, context) - faults = self.parameterAsSource(parameters, self.INPUT_FAULTS, context) - strati_column = self.parameterAsSource(parameters, self.INPUT_STRATI_COLUMN, context) - - geology = qgsLayerToGeoDataFrame(geology) - faults = qgsLayerToGeoDataFrame(faults) if faults else None - - feedback.pushInfo("Extracting Basal Contacts...") - contact_extractor = ContactExtractor(geology, faults, feedback) - contact_extractor.extract_basal_contacts(strati_column) - - basal_contacts = GeoDataFrameToQgsLayer( - self, - contact_extractor.basal_contacts, - parameters=parameters, - context=context, - feedback=feedback, - ) - return {self.OUTPUT: basal_contacts} - - def createInstance(self) -> QgsProcessingAlgorithm: - """Create a new instance of the algorithm.""" - return self.__class__() # BasalContactsAlgorithm() - - - -class StructuralPointAlgorithm(QgsProcessingAlgorithm): - """Processing algorithm for thickness calculations.""" - - - INPUT_UNITS = 'UNITS' - INPUT_STRATI_COLUMN = 'STRATIGRAPHIC_COLUMN' - INPUT_BASAL_CONTACTS = 'BASAL_CONTACTS' - INPUT_STRUCTURE_DATA = 'STRUCTURE_DATA' - INPUT_GEOLOGY = 'GEOLOGY' - INPUT_SAMPLED_CONTACTS = 'SAMPLED_CONTACTS' - - OUTPUT = "THICKNESS" - - def name(self) -> str: - """Return the algorithm name.""" - return "thickness_calculator" - - def displayName(self) -> str: - """Return the algorithm display name.""" - return "Loop3d: Thickness Calculator" - - def group(self) -> str: - """Return the algorithm group name.""" - return "Loop3d" - - def groupId(self) -> str: - """Return the algorithm group ID.""" - return "loop3d" - - def initAlgorithm(self, config: Optional[dict[str, Any]] = None) -> None: - """Initialize the algorithm parameters.""" - - self.addParameter( - QgsProcessingParameterFeatureSource( - self.INPUT_GEOLOGY, - "GEOLOGY", - [QgsProcessing.TypeVectorPolygon], + QgsProcessingParameterMatrix( + name=self.INPUT_STRATI_COLUMN, + description="Stratigraphic Order", + headers=["Unit"], + numberRows=0, + defaultValue=[] ) ) self.addParameter( QgsProcessingParameterFeatureSource( - self.INPUT_FAULTS, - "FAULTS", - [QgsProcessing.TypeVectorLine], - optional=True, + self.INPUT_SAMPLED_CONTACTS, + "SAMPLED_CONTACTS", + [QgsProcessing.TypeVectorPoint], ) ) - self.addParameter( QgsProcessingParameterFeatureSource( - self.INPUT_STRATI_COLUMN, - "STRATIGRAPHIC_COLUMN", - [QgsProcessing.TypeVectorLine], + self.INPUT_STRUCTURE_DATA, + "STRUCTURE_DATA", + [QgsProcessing.TypeVectorPoint], ) ) - self.addParameter( QgsProcessingParameterFeatureSink( self.OUTPUT, - "Basal Contacts", + "Thickness", ) ) @@ -319,25 +156,66 @@ def processAlgorithm( feedback: QgsProcessingFeedback, ) -> dict[str, Any]: - geology = self.parameterAsSource(parameters, self.INPUT_GEOLOGY, context) - faults = self.parameterAsSource(parameters, self.INPUT_FAULTS, context) - strati_column = self.parameterAsSource(parameters, self.INPUT_STRATI_COLUMN, context) - - geology = qgsLayerToGeoDataFrame(geology) - faults = qgsLayerToGeoDataFrame(faults) if faults else None - - feedback.pushInfo("Extracting Basal Contacts...") - contact_extractor = ContactExtractor(geology, faults, feedback) - contact_extractor.extract_basal_contacts(strati_column) - - basal_contacts = GeoDataFrameToQgsLayer( + feedback.pushInfo("Initialising Thickness Calculation Algorithm...") + thickness_type = self.parameterAsEnum(parameters, self.INPUT_THICKNESS_CALCULATOR_TYPE, context) + dtm_data = self.parameterAsSource(parameters, self.INPUT_DTM, context) + bounding_box = self.parameterAsEnum(parameters, self.INPUT_BOUNDING_BOX, context) + max_line_length = self.parameterAsNumber(parameters, self.INPUT_MAX_LINE_LENGTH, context) + units = self.parameterAsSource(parameters, self.INPUT_UNITS, context) + basal_contacts = self.parameterAsSource(parameters, self.INPUT_BASAL_CONTACTS, context) + geology_data = self.parameterAsSource(parameters, self.INPUT_GEOLOGY, context) + stratigraphic_order = self.parameterAsMatrix(parameters, self.INPUT_STRATI_COLUMN, context) + structure_data = self.parameterAsSource(parameters, self.INPUT_STRUCTURE_DATA, context) + sampled_contacts = self.parameterAsSource(parameters, self.INPUT_SAMPLED_CONTACTS, context) + + # convert layers to dataframe or geodataframe + geology_data = qgsLayerToGeoDataFrame(geology_data) + units = qgsLayerToDataFrame(units) + basal_contacts = qgsLayerToGeoDataFrame(basal_contacts) + structure_data = qgsLayerToDataFrame(structure_data) + sampled_contacts = qgsLayerToDataFrame(sampled_contacts) + + feedback.pushInfo("Calculating unit thicknesses...") + + if thickness_type == "InterpolatedStructure": + thickness_calculator = InterpolatedStructure( + dtm_data=dtm_data, + bounding_box=bounding_box, + ) + thickness_calculator.compute( + units, + stratigraphic_order, + basal_contacts, + structure_data, + geology_data, + sampled_contacts + ) + + if thickness_type == "StructuralPoint": + thickness_calculator = StructuralPoint( + dtm_data=dtm_data, + bounding_box=bounding_box, + max_line_length=max_line_length, + ) + thickness_calculator.compute( + units, + stratigraphic_order, + basal_contacts, + structure_data, + geology_data, + sampled_contacts + ) + + #TODO: convert thicknesses dataframe to qgs layer + thicknesses = dataframeToQgsLayer( self, - contact_extractor.basal_contacts, + # contact_extractor.basal_contacts, parameters=parameters, context=context, feedback=feedback, ) - return {self.OUTPUT: basal_contacts} + + return {self.OUTPUT: thicknesses[1]} def createInstance(self) -> QgsProcessingAlgorithm: """Create a new instance of the algorithm."""