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 }} diff --git a/map2loop/main/vectorLayerWrapper.py b/map2loop/main/vectorLayerWrapper.py new file mode 100644 index 0000000..d68ed3d --- /dev/null +++ b/map2loop/main/vectorLayerWrapper.py @@ -0,0 +1,451 @@ +# PyQGIS / PyQt imports + +from qgis.core import ( + QgsRaster, + QgsFields, + QgsField, + QgsFeature, + QgsGeometry, + QgsWkbTypes, + QgsCoordinateReferenceSystem, + QgsFeatureSink, + QgsProcessingException + ) + +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 + + + +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) + +def GeoDataFrameToQgsLayer(qgs_algorithm, geodataframe, parameters, context, output_key, feedback=None): + """ + 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 } + """ + + 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 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 = QgsCoordinateReferenceSystem.fromWkt(geodataframe.crs.to_wkt()) + except Exception: + try: + epsg = geodataframe.crs.to_epsg() + if epsg: + crs = QgsCoordinateReferenceSystem.fromEpsgId(int(epsg)) + except Exception: + pass + + # --- build QGIS fields from pandas dtypes + fields = QgsFields() + 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): + 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: + 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 + attrs = [] + for col in non_geom_cols: + val = row[col] + if isinstance(val, np.generic): + try: + val = val.item() + except Exception: + pass + if pd.api.types.is_datetime64_any_dtype(geodataframe[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)) + + sink.addFeature(f, QgsFeatureSink.FastInsert) + + if total: + feedback.setProgress(int(100.0 * (i + 1) / total)) + + 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 dataframeToQgsLayer( + 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 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..534958e --- /dev/null +++ b/map2loop/processing/algorithms/extract_basal_contacts.py @@ -0,0 +1,119 @@ +""" +*************************************************************************** +* * +* 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 import ContactExtractor + + +class BasalContactsAlgorithm(QgsProcessingAlgorithm): + """Processing algorithm to create basal contacts.""" + + + INPUT_GEOLOGY = 'GEOLOGY' + INPUT_FAULTS = 'FAULTS' + INPUT_STRATI_COLUMN = 'STRATIGRAPHIC_COLUMN' + OUTPUT = "BASAL_CONTACTS" + + def name(self) -> str: + """Return the algorithm name.""" + return "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, + "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() diff --git a/map2loop/processing/algorithms/sampler.py b/map2loop/processing/algorithms/sampler.py new file mode 100644 index 0000000..3feb3e6 --- /dev/null +++ b/map2loop/processing/algorithms/sampler.py @@ -0,0 +1,191 @@ +""" +*************************************************************************** +* * +* 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 +from qgis.PyQt.QtCore import QMetaType + +# QGIS imports +from qgis.core import ( + QgsFeatureSink, + QgsProcessing, + QgsProcessingAlgorithm, + QgsProcessingContext, + QgsProcessingException, + QgsProcessingFeedback, + QgsProcessingParameterFeatureSink, + QgsProcessingParameterFeatureSource, + QgsProcessingParameterString, + QgsProcessingParameterNumber, + QgsField, + QgsFeature, + QgsGeometry, + QgsPointXY, + QgsVectorLayer +) +# Internal imports +from ...main.vectorLayerWrapper import qgsLayerToGeoDataFrame +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 == "decimator": + feedback.pushInfo("Sampling...") + sampler = SamplerDecimator(decimation=decimation, dtm_data=dtm, geology_data=geology, feedback=feedback) + samples = sampler.sample(spatial_data) + + if sampler_type == "spacing": + feedback.pushInfo("Sampling...") + sampler = SamplerSpacing(spacing=spacing, dtm_data=dtm, geology_data=geology, feedback=feedback) + samples = sampler.sample(spatial_data) + + + # 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.""" + return self.__class__() # SamplerAlgorithm() \ No newline at end of file diff --git a/map2loop/processing/algorithms/sorter.py b/map2loop/processing/algorithms/sorter.py new file mode 100644 index 0000000..a159855 --- /dev/null +++ b/map2loop/processing/algorithms/sorter.py @@ -0,0 +1,205 @@ +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, + QgsWkbTypes +) + +# ──────────────────────────────────────────────── +# 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 diff --git a/map2loop/processing/algorithms/thickness_calculator.py b/map2loop/processing/algorithms/thickness_calculator.py new file mode 100644 index 0000000..f56b09c --- /dev/null +++ b/map2loop/processing/algorithms/thickness_calculator.py @@ -0,0 +1,222 @@ +""" +*************************************************************************** +* * +* 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, + QgsProcessingParameterEnum, + QgsProcessingParameterNumber, + QgsProcessingParameterField, + QgsProcessingParameterMatrix +) +# Internal imports +from ...main.vectorLayerWrapper import qgsLayerToGeoDataFrame, GeoDataFrameToQgsLayer, qgsLayerToDataFrame, dataframeToQgsLayer +from map2loop.map2loop.thickness_calculator import InterpolatedStructure, StructuralPoint + + +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' + 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( + QgsProcessingParameterEnum( + self.INPUT_THICKNESS_CALCULATOR_TYPE, + "Thickness Calculator Type", + options=['InterpolatedStructure','StructuralPoint'], + allowMultiple=False, + ) + ) + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_DTM, + "DTM", + [QgsProcessing.TypeVectorRaster], + ) + ) + self.addParameter( + 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_UNITS, + "Units", + [QgsProcessing.TypeVectorLine], + ) + ) + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_BASAL_CONTACTS, + "Basal Contacts", + [QgsProcessing.TypeVectorLine], + ) + ) + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_GEOLOGY, + "GEOLOGY", + [QgsProcessing.TypeVectorPolygon], + ) + ) + self.addParameter( + QgsProcessingParameterMatrix( + name=self.INPUT_STRATI_COLUMN, + description="Stratigraphic Order", + headers=["Unit"], + numberRows=0, + defaultValue=[] + ) + ) + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_SAMPLED_CONTACTS, + "SAMPLED_CONTACTS", + [QgsProcessing.TypeVectorPoint], + ) + ) + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_STRUCTURE_DATA, + "STRUCTURE_DATA", + [QgsProcessing.TypeVectorPoint], + ) + ) + self.addParameter( + QgsProcessingParameterFeatureSink( + self.OUTPUT, + "Thickness", + ) + ) + + def processAlgorithm( + self, + parameters: dict[str, Any], + context: QgsProcessingContext, + feedback: QgsProcessingFeedback, + ) -> dict[str, Any]: + + 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, + parameters=parameters, + context=context, + feedback=feedback, + ) + + return {self.OUTPUT: thicknesses[1]} + + def createInstance(self) -> QgsProcessingAlgorithm: + """Create a new instance of the algorithm.""" + return self.__class__() # BasalContactsAlgorithm() \ No newline at end of file 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: