diff --git a/.github/labeler.yml b/.github/labeler.yml index a16bed4..bf6729d 100644 --- a/.github/labeler.yml +++ b/.github/labeler.yml @@ -45,7 +45,6 @@ tooling: - .pre-commit-config.yaml - setup.cfg - UI: - head-branch: - ^ui diff --git a/.github/workflows/linter.yml b/.github/workflows/linter.yml index ef30cb0..78979a2 100644 --- a/.github/workflows/linter.yml +++ b/.github/workflows/linter.yml @@ -19,7 +19,6 @@ env: permissions: contents: write - jobs: lint-py: name: Python 🐍 diff --git a/.github/workflows/package_and_release.yml b/.github/workflows/package_and_release.yml new file mode 100644 index 0000000..63632c3 --- /dev/null +++ b/.github/workflows/package_and_release.yml @@ -0,0 +1,168 @@ +name: "πŸ“¦ Package & πŸš€ Release" + +env: + PROJECT_FOLDER: plugin_map2loop + PYTHON_VERSION: 3.9 + +on: + push: + branches: + - main + paths: + - .github/workflows/package_and_release.yml + - 'docs/**/*' + - "plugin_map2loop/**/*.py" + - "plugin_map2loop/metadata.txt" + tags: + - "*" + +# Allow one concurrent deployment per branch/pr +concurrency: + group: ${{ github.workflow }}-${{ github.event.pull_request.number || github.ref }} + cancel-in-progress: true + +jobs: + translation: + name: "πŸ’¬ i18n compilation" + runs-on: ubuntu-latest + + steps: + - name: Get source code + uses: actions/checkout@v4 + + - name: Setup Python + uses: actions/setup-python@v5 + with: + python-version: ${{ env.PYTHON_VERSION }} + + - name: Install translation requirements + run: | + sudo apt update + sudo apt install qt5-qmake qttools5-dev-tools + python3 -m pip install -U pyqt5-tools + + - name: Update translations + run: | + python3 scripts/generate_translation_profile.py + pylupdate5 -noobsolete -verbose ${{ env.PROJECT_FOLDER }}/resources/i18n/plugin_translation.pro + + - name: Compile translations + run: lrelease ${{ env.PROJECT_FOLDER }}/resources/i18n/*.ts + + - uses: actions/upload-artifact@v4 + with: + name: translations-build + path: ${{ env.PROJECT_FOLDER }}/**/*.qm + if-no-files-found: error + + # -- NO TAGS ---------------------------------------------------------------------- + packaging: + name: "πŸ“¦ Packaging plugin" + runs-on: ubuntu-latest + needs: + - translation + + if: ${{ !startsWith(github.ref, 'refs/tags/') }} + + steps: + - name: Checkout + uses: actions/checkout@v4 + + - name: Setup Python + uses: actions/setup-python@v5 + with: + cache: "pip" + cache-dependency-path: "requirements/packaging.txt" + python-version: ${{ env.PYTHON_VERSION }} + + - name: Install dependencies + run: | + python -m pip install -U pip setuptools wheel + python -m pip install -U -r requirements/packaging.txt + + - name: Download translations + uses: actions/download-artifact@v4 + with: + name: translations-build + path: ${{ env.PROJECT_FOLDER }} + + - name: Amend gitignore to include compiled translations and add it to tracked files + run: | + # include compiled translations + sed -i "s|^*.qm.*| |" .gitignore + + # git add full project + git add ${{ env.PROJECT_FOLDER }}/ + + - name: Package the latest version + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + run: | + qgis-plugin-ci package latest \ + --allow-uncommitted-changes \ + --plugin-repo-url $(gh api "repos/$GITHUB_REPOSITORY/pages" --jq '.html_url') + + - uses: actions/upload-artifact@v4 + with: + name: ${{ env.PROJECT_FOLDER }}-latest + path: | + plugins.xml + ${{ env.PROJECT_FOLDER }}.*.zip + if-no-files-found: error + + # -- ONLY TAGS ---------------------------------------------------------------------- + release: + name: "πŸš€ Release on tag" + runs-on: ubuntu-latest + permissions: + contents: write + needs: + - translation + + if: startsWith(github.ref, 'refs/tags/') + + steps: + - name: Checkout + uses: actions/checkout@v4 + + - name: Setup Python + uses: actions/setup-python@v5 + with: + cache: "pip" + cache-dependency-path: "requirements/packaging.txt" + python-version: ${{ env.PYTHON_VERSION }} + + - name: Install project requirements + run: | + python -m pip install -U pip setuptools wheel + python -m pip install -U -r requirements/packaging.txt + + - name: Download translations + uses: actions/download-artifact@v4 + with: + name: translations-build + path: ${{ env.PROJECT_FOLDER }} + + - name: Amend gitignore to include compiled translations and add it to tracked files + run: | + # include compiled translations + sed -i "s|^*.qm.*| |" .gitignore + + # git add full project + git add ${{ env.PROJECT_FOLDER }}/ + + - name: Create GitHub Release + uses: softprops/action-gh-release@v2 + with: + fail_on_unmatched_files: true + generate_release_notes: true + + - name: Deploy plugin + run: >- + qgis-plugin-ci + release ${GITHUB_REF/refs\/tags\//} + --allow-uncommitted-changes + --create-plugin-repo + --github-token ${{ secrets.GITHUB_TOKEN }} + --osgeo-username ${{ secrets.OSGEO_USER }} + --osgeo-password ${{ secrets.OSGEO_PASSWORD }} \ No newline at end of file diff --git a/.github/workflows/release-please.yml b/.github/workflows/release-please.yml index 5805954..8862ca3 100644 --- a/.github/workflows/release-please.yml +++ b/.github/workflows/release-please.yml @@ -1,4 +1,3 @@ - on: push: branches: diff --git a/.release-please-manifest.json b/.release-please-manifest.json index d28b6e6..3f6ead5 100644 --- a/.release-please-manifest.json +++ b/.release-please-manifest.json @@ -1,4 +1,4 @@ { ".": "0.0.1", "loopstructural": "0.1.11" -} +} \ No newline at end of file diff --git a/AGENTS.md b/AGENTS.md new file mode 100644 index 0000000..58b90f9 --- /dev/null +++ b/AGENTS.md @@ -0,0 +1,93 @@ +# AGENTS.md + +## Purpose + +This document outlines the architectural and development principles for contributors to the `plugin_loopstructural` QGIS plugin. The plugin provides a thin, modular interface to the `map2loop` and `LoopStructural` libraries, enabling geological modeling workflows within QGIS. + +--- + +## Design Philosophy + +### πŸ”Ή Thin Interface Layer +- The plugin **must not** reimplement or duplicate functionality from `map2loop` or `LoopStructural`. +- All core logic and enhancements should be contributed upstream to the respective libraries. +- The plugin should focus on **UI integration**, **data flow orchestration**, and **user interaction**. + +### πŸ”Ή Modularity +- UI components (dialogs, panels, actions) live under the `loopstructural/gui/` package and should be encapsulated in their own classes. +- Business logic and orchestration are located in `loopstructural/main/` and `loopstructural/toolbelt/` where adapters and services wrap external libraries. +- Processing algorithms and QGIS provider integration are in `loopstructural/processing/` and should be isolated from UI code. +- Avoid tight coupling between components. Use signals/slots or event-driven patterns where appropriate. + +### πŸ”Ή Object-Oriented Design +- Use classes with clear responsibilities and interfaces. +- Prefer composition over inheritance unless subclassing is semantically appropriate. +- Encapsulate interactions with external libraries in dedicated adapter or service classes (e.g., `loopstructural.main.Map2LoopService`, `loopstructural.main.LoopStructuralRunner`). + +--- + +## Development Guidelines + +### βœ… Code Quality +- All code must pass the repository's pre-commit checks (formatting, linting, import sorting). +- Use type hints and docstrings for all public methods and classes. +- Follow PEP8 and QGIS plugin development best practices. + +### πŸ§ͺ Testing +- All new code must include **unit tests** and, where applicable, **integration tests**. +- Tests live under the `tests/` package and are runnable with `pytest`. +- Mock external dependencies (`map2loop`, `LoopStructural`) in unit tests. + +### 🧩 Current Plugin Structure + +``` +plugin_loopstructural/ +β”œβ”€β”€ loopstructural/ # plugin package +β”‚ β”œβ”€β”€ __init__.py +β”‚ β”œβ”€β”€ __about__.py +β”‚ β”œβ”€β”€ plugin_main.py # QGIS plugin entry and bootstrap +β”‚ β”œβ”€β”€ gui/ # UI dialogs, widgets, and panels +β”‚ β”œβ”€β”€ main/ # controllers, managers, adapters (service layer) +β”‚ β”œβ”€β”€ processing/ # QGIS processing provider and algorithms +β”‚ β”œβ”€β”€ toolbelt/ # utilities, env parsing, preferences, logging +β”‚ β”œβ”€β”€ resources/ # icons, translations, help files +β”‚ └── requirements.txt +β”œβ”€β”€ docs/ +β”œβ”€β”€ requirements/ +β”œβ”€β”€ tests/ +└── README.md +``` + +Notes on mapping older concepts: +- What used to be called `services/` and `controllers/` is implemented across `loopstructural/main/` and `loopstructural/toolbelt/`. +- UI remains in `loopstructural/gui/` (dialogs, `.ui` files, widget classes). +- Processing-specific code and QGIS provider live under `loopstructural/processing/`. + +--- + +## Contribution Workflow + +1. Fork the repository and create a feature branch. +2. Implement changes following the design and code quality guidelines. +3. Add or update tests under `tests/` and ensure they run with `pytest`. +4. Run pre-commit hooks (e.g. `pre-commit run --all-files`) and ensure all checks pass. +5. Submit a pull request with a clear description of the changes and rationale. Link to upstream libraries if behavior is moved upstream. + +--- + +## Future Enhancements + +- Support for asynchronous or background processing of long-running tasks (consider using QGIS background task framework). +- Improved error handling and user feedback. +- Internationalization (i18n) support and keeping `.ts`/.qm translation files in `loopstructural/resources/i18n/`. +- Better separation of concerns between UI, processing algorithms and adapters to facilitate unit testing. + +--- + +## Contact + +For questions or contributions to the upstream libraries: +- `map2loop` +- `LoopStructural` + +--- diff --git a/LICENSE b/LICENSE index c2811fc..5f5da4d 100644 --- a/LICENSE +++ b/LICENSE @@ -2,7 +2,7 @@ GNU GENERAL PUBLIC LICENSE Version 2, June 1991 - Copyright (c) 2024, Lachlan GROSE / QGIS + Copyright (c) 2025, Lachlan GROSE / QGIS Copyright (C) 1989, 1991 Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA Everyone is permitted to copy and distribute verbatim copies diff --git a/docs/conf.py b/docs/conf.py index f652ebb..86dc9ca 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -1,7 +1,6 @@ #!python3 -"""Configuration for project documentation using Sphinx. -""" +"""Configuration for project documentation using Sphinx.""" # standard import sys diff --git a/docs/development/environment.md b/docs/development/environment.md index f8f0ee2..f63afc4 100644 --- a/docs/development/environment.md +++ b/docs/development/environment.md @@ -8,7 +8,6 @@ Typically on Ubuntu: # create virtual environment linking to system packages (for pyqgis) python3 -m venv .venv --system-site-packages source .venv/bin/activate - # bump dependencies inside venv python -m pip install -U pip python -m pip install -U -r requirements/development.txt diff --git a/docs/static/dev_qgis_enable_plugin.png b/docs/static/dev_qgis_enable_plugin.png new file mode 100644 index 0000000..ec9cbc3 Binary files /dev/null and b/docs/static/dev_qgis_enable_plugin.png differ diff --git a/docs/static/dev_qgis_set_pluginpath_envvar.png b/docs/static/dev_qgis_set_pluginpath_envvar.png new file mode 100644 index 0000000..9c0ecc1 Binary files /dev/null and b/docs/static/dev_qgis_set_pluginpath_envvar.png differ diff --git a/loopstructural/main/vectorLayerWrapper.py b/loopstructural/main/vectorLayerWrapper.py index f84f6d0..4ec74d8 100644 --- a/loopstructural/main/vectorLayerWrapper.py +++ b/loopstructural/main/vectorLayerWrapper.py @@ -1,6 +1,92 @@ -import pandas as pd +import os +import tempfile + import geopandas as gpd -from qgis.core import QgsRaster, QgsWkbTypes +import numpy as np +import pandas as pd + +# PyQGIS / PyQt imports +from osgeo import gdal +from qgis import processing +from qgis.core import ( + QgsCoordinateReferenceSystem, + QgsCoordinateTransform, + QgsFeature, + QgsFeatureSink, + QgsField, + QgsFields, + QgsGeometry, + QgsPoint, + QgsPointXY, + QgsProcessingException, + QgsProject, + QgsRaster, + QgsRasterLayer, + QgsWkbTypes, +) +from qgis.PyQt.QtCore import QDateTime, QVariant +from shapely.geometry import LineString, MultiLineString, MultiPoint, MultiPolygon, Point, Polygon +from shapely.wkb import loads as wkb_loads + + +def qgsRasterToGdalDataset(rlayer: QgsRasterLayer): + """ + Convert a QgsRasterLayer to an osgeo.gdal.Dataset (read-only). + If the raster is non-file-based (e.g. WMS/WCS/virtual), we create a temp GeoTIFF via gdal:translate. + Returns a gdal.Dataset or None. + """ + if rlayer is None or not rlayer.isValid(): + return None + + # Try direct open on file-backed layers + candidates = [] + try: + candidates.append(rlayer.source()) + except Exception: + pass + try: + if rlayer.dataProvider(): + candidates.append(rlayer.dataProvider().dataSourceUri()) + except Exception: + pass + + tried = set() + for uri in candidates: + if not uri: + continue + if uri in tried: + continue + tried.add(uri) + + # Strip QGIS pipe options: "path.tif|layername=..." β†’ "path.tif" + base_uri = uri.split("|")[0] + + # Some providers store β€œSUBDATASET:” URIs; gdal.OpenEx can usually handle them directly. + ds = gdal.OpenEx(base_uri, gdal.OF_RASTER | gdal.OF_READONLY) + if ds is not None: + return ds + + # If we’re here, it’s likely non-file-backed. Export to a temp GeoTIFF. + tmpdir = tempfile.gettempdir() + tmp_path = os.path.join(tmpdir, f"m2l_dtm_{rlayer.id()}.tif") + + # Use GDAL Translate via QGIS processing (avoids CRS pitfalls) + processing.run( + "gdal:translate", + { + "INPUT": rlayer, # QGIS accepts the layer object here + "TARGET_CRS": None, + "NODATA": None, + "COPY_SUBDATASETS": False, + "OPTIONS": "", + "EXTRA": "", + "DATA_TYPE": 0, # Use input data type + "OUTPUT": tmp_path, + }, + ) + + ds = gdal.OpenEx(tmp_path, gdal.OF_RASTER | gdal.OF_READONLY) + return ds def qgsLayerToGeoDataFrame(layer) -> gpd.GeoDataFrame: @@ -17,69 +103,660 @@ def qgsLayerToGeoDataFrame(layer) -> gpd.GeoDataFrame: continue data['geometry'].append(geom) for f in fields: - data[f.name()].append(feature[f.name()]) - return gpd.GeoDataFrame(data, crs=layer.crs().authid()) + if f.type() == QVariant.String: + data[f.name()].append(str(feature[f.name()])) + else: + data[f.name()].append(feature[f.name()]) + return gpd.GeoDataFrame(data, crs=layer.sourceCrs().authid()) + + +def qgsLayerToDataFrame(src, dtm=None) -> pd.DataFrame: + """ + Convert a vector layer or processing feature source to a pandas DataFrame. + Samples geometry using points or vertices of lines/polygons. + Optionally samples Z from a DTM raster. + + :param src: QgsVectorLayer or QgsProcessingFeatureSource + :param dtm: QgsRasterLayer or None + :return: pd.DataFrame with columns: X, Y, Z, and all layer fields + """ + + if src is None: + return None + + # --- Resolve fields and source CRS (works for both layer and feature source) --- + fields = src.fields() if hasattr(src, "fields") else None + if fields is None: + # Fallback: take fields from first feature if needed + feat_iter = src.getFeatures() + try: + first = next(feat_iter) + except StopIteration: + return pd.DataFrame(columns=["X", "Y", "Z"]) + fields = first.fields() + # Rewind iterator by building a new one + feats = [first] + list(src.getFeatures()) + else: + feats = src.getFeatures() + + # Get source CRS + if hasattr(src, "crs"): + src_crs = src.crs() + elif hasattr(src, "sourceCrs"): + src_crs = src.sourceCrs() + else: + src_crs = None + + # --- Prepare optional transform to DTM CRS for sampling --- + to_dtm = None + if dtm is not None and src_crs is not None and dtm.crs().isValid() and src_crs.isValid(): + if src_crs != dtm.crs(): + to_dtm = QgsCoordinateTransform(src_crs, dtm.crs(), QgsProject.instance()) + + # --- Helper: sample Z from DTM (returns float or -9999) --- + def sample_dtm_xy(x, y): + if dtm is None: + return 0.0 + # Transform coordinate if needed + if to_dtm is not None: + try: + from qgis.core import QgsPointXY + + x, y = to_dtm.transform(QgsPointXY(x, y)) + except Exception: + return -9999.0 + from qgis.core import QgsPointXY + + ident = dtm.dataProvider().identify(QgsPointXY(x, y), QgsRaster.IdentifyFormatValue) + if not ident.isValid(): + return -9999.0 + res = ident.results() + if not res: + return -9999.0 + # take first band value (band keys are 1-based) + try: + # Prefer band 1 if present + return float(res.get(1, next(iter(res.values())))) + except Exception: + return -9999.0 + + # --- Geometry -> list of vertices (QgsPoint or QgsPointXY) --- + def vertices_from_geometry(geom): + if geom is None or geom.isEmpty(): + return [] + gtype = QgsWkbTypes.geometryType(geom.wkbType()) + is_multi = QgsWkbTypes.isMultiType(geom.wkbType()) + + if gtype == QgsWkbTypes.PointGeometry: + if is_multi: + return list(geom.asMultiPoint()) + else: + return [geom.asPoint()] + + elif gtype == QgsWkbTypes.LineGeometry: + pts = [] + if is_multi: + for line in geom.asMultiPolyline(): + pts.extend(line) + else: + pts.extend(geom.asPolyline()) + return pts + + elif gtype == QgsWkbTypes.PolygonGeometry: + pts = [] + if is_multi: + mpoly = geom.asMultiPolygon() + for poly in mpoly: + for ring in poly: # exterior + interior rings + pts.extend(ring) + else: + poly = geom.asPolygon() + for ring in poly: + pts.extend(ring) + return pts + + # Other geometry types not handled + return [] + # --- Build rows safely (one dict per sampled point) --- + rows = [] + field_names = [f.name() for f in fields] -def qgsLayerToDataFrame(layer, dtm) -> pd.DataFrame: - """Convert a vector layer to a pandas DataFrame by sampling geometries. + for f in feats: + geom = f.geometry() + pts = vertices_from_geometry(geom) - Samples point geometries or the vertices of line geometries and optionally - queries a DTM for Z values. + if not pts: + # If you want to keep attribute rows even when no vertices: uncomment below + # row = {name: f[name] for name in field_names} + # row.update({"X": None, "Y": None, "Z": None}) + # rows.append(row) + continue + + # Cache attributes once per feature and reuse for each sampled point + base_attrs = {name: f[name] for name in field_names} + + for p in pts: + # QgsPoint vs QgsPointXY both have x()/y() + x, y = float(p.x()), float(p.y()) + z = sample_dtm_xy(x, y) + + row = {"X": x, "Y": y, "Z": z} + row.update(base_attrs) + rows.append(row) + + # Create DataFrame; if empty, return with expected columns + if not rows: + cols = ["X", "Y", "Z"] + field_names + return pd.DataFrame(columns=cols) + + return pd.DataFrame.from_records(rows) + + +def GeoDataFrameToQgsLayer( + qgs_algorithm, geodataframe, parameters, context, output_key, feedback=None +): + """ + Write a GeoPandas GeoDataFrame directly to a QGIS Processing FeatureSink. Parameters ---------- - layer : QgsVectorLayer - Input QGIS vector layer to sample. - dtm : QgsRaster or None - Digital Terrain Model used to evaluate Z values for points (optional). + alg : QgsProcessingAlgorithm (self) + gdf : geopandas.GeoDataFrame + parameters : dict (from processAlgorithm) + context : QgsProcessingContext + output_key : str (e.g. self.OUTPUT) + feedback : QgsProcessingFeedback | None Returns ------- - pd.DataFrame - DataFrame with columns `X`, `Y`, `Z` and one column per layer field. + str : dest_id to return from processAlgorithm, e.g. { output_key: dest_id } """ - 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] + + 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: - 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) + 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 + + +def matrixToDict(matrix, headers=("minx", "miny", "maxx", "maxy")) -> dict: + """ + Convert a QgsProcessingParameterMatrix value to a dict with float values. + Accepts: [[minx,miny,maxx,maxy]] or [minx,miny,maxx,maxy]. + Raises a clear error if an enum index (int) was passed by mistake. + """ + # Guard: common mistake β†’ using parameterAsEnum + if isinstance(matrix, int): + raise QgsProcessingException( + "Bounding Box was read with parameterAsEnum (got an int). " + "Use parameterAsMatrix for QgsProcessingParameterMatrix." + ) + + if matrix is None: + raise QgsProcessingException("Bounding box matrix is None.") + + # Allow empty string from settings/defaults + if isinstance(matrix, str) and not matrix.strip(): + raise QgsProcessingException("Bounding box matrix is empty.") + + # Accept single-row matrix or flat list + if isinstance(matrix, (list, tuple)): + if matrix and isinstance(matrix[0], (list, tuple)): + row = matrix[0] + else: + row = matrix + else: + # last resort: try comma-separated string "minx,miny,maxx,maxy" + if isinstance(matrix, str) and "," in matrix: + row = [v.strip() for v in matrix.split(",")] + else: + raise QgsProcessingException(f"Unrecognized bounding box value: {type(matrix)}") + + if len(row) < 4: + raise QgsProcessingException(f"Bounding box needs 4 numbers, got {len(row)}: {row}") + + def _to_float(v): + if isinstance(v, str): + v = v.strip() + return float(v) + + vals = list(map(_to_float, row[:4])) + bbox = dict(zip(headers, vals)) + + if not (bbox["minx"] < bbox["maxx"] and bbox["miny"] < bbox["maxy"]): + raise QgsProcessingException( + f"Invalid bounding box: {bbox} (expect minx (value, ok) + d, ok = val.toDouble() + return float(d) if ok else None + # native int/float + if isinstance(val, (int, float)): + return float(val) + # fallback conversion attempt + try: + return float(val) + except Exception: + return None diff --git a/loopstructural/plugin_main.py b/loopstructural/plugin_main.py index 795eef1..19d01c8 100644 --- a/loopstructural/plugin_main.py +++ b/loopstructural/plugin_main.py @@ -35,6 +35,9 @@ from loopstructural.gui.loop_widget import LoopWidget from loopstructural.main.data_manager import ModellingDataManager from loopstructural.main.model_manager import GeologicalModelManager +from loopstructural.processing import ( + Map2LoopProvider, +) from loopstructural.toolbelt import PlgLogger, PlgOptionsManager # ############################################################################ @@ -137,6 +140,7 @@ def initGui(self): # -- Menu self.iface.addPluginToMenu(__title__, self.action_settings) self.iface.addPluginToMenu(__title__, self.action_help) + self.initProcessing() # -- Help menu @@ -266,6 +270,11 @@ def tr(self, message: str) -> str: """ return QCoreApplication.translate(self.__class__.__name__, message) + def initProcessing(self): + """Initialize the processing provider.""" + self.provider = Map2LoopProvider() + QgsApplication.processingRegistry().addProvider(self.provider) + def unload(self): """Clean up when plugin is disabled or uninstalled.""" # -- Clean up dock widgets @@ -285,6 +294,8 @@ def unload(self): # self.iface.removeMenu(self.menu) # -- Clean up preferences panel in QGIS settings self.iface.unregisterOptionsWidgetFactory(self.options_factory) + # -- Unregister processing + QgsApplication.processingRegistry().removeProvider(self.provider) # remove from QGIS help/extensions menu if self.action_help_plugin_menu_documentation: diff --git a/loopstructural/processing/__init__.py b/loopstructural/processing/__init__.py index 6b649c5..c335d6f 100644 --- a/loopstructural/processing/__init__.py +++ b/loopstructural/processing/__init__.py @@ -4,4 +4,4 @@ Contains provider definitions and registration utilities for the QGIS processing framework. """ -from .provider import LoopstructuralProvider +from .provider import LoopstructuralProvider, Map2LoopProvider diff --git a/loopstructural/processing/algorithms/__init__.py b/loopstructural/processing/algorithms/__init__.py new file mode 100644 index 0000000..08d76a0 --- /dev/null +++ b/loopstructural/processing/algorithms/__init__.py @@ -0,0 +1,5 @@ +from .extract_basal_contacts import BasalContactsAlgorithm +from .sorter import StratigraphySorterAlgorithm +from .user_defined_sorter import UserDefinedStratigraphyAlgorithm +from .thickness_calculator import ThicknessCalculatorAlgorithm +from .sampler import SamplerAlgorithm diff --git a/loopstructural/processing/algorithms/extract_basal_contacts.py b/loopstructural/processing/algorithms/extract_basal_contacts.py new file mode 100644 index 0000000..6446b13 --- /dev/null +++ b/loopstructural/processing/algorithms/extract_basal_contacts.py @@ -0,0 +1,194 @@ +""" +*************************************************************************** +* * +* 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, + QgsProcessingParameterMapLayer, + QgsProcessingParameterString, + QgsProcessingParameterField, + QgsProcessingParameterMatrix, + QgsVectorLayer, + QgsSettings +) +# Internal imports +from ...main.vectorLayerWrapper import qgsLayerToGeoDataFrame, GeoDataFrameToQgsLayer +from map2loop.contact_extractor import ContactExtractor + + +class BasalContactsAlgorithm(QgsProcessingAlgorithm): + """Processing algorithm to create basal contacts.""" + + + INPUT_GEOLOGY = 'GEOLOGY' + INPUT_FAULTS = 'FAULTS' + INPUT_STRATI_COLUMN = 'STRATIGRAPHIC_COLUMN' + INPUT_IGNORE_UNITS = 'IGNORE_UNITS' + OUTPUT = "BASAL_CONTACTS" + ALL_CONTACTS = "ALL_CONTACTS" + + def name(self) -> str: + """Return the algorithm name.""" + return "basal_contacts" + + def displayName(self) -> str: + """Return the algorithm display name.""" + return "Basal Contacts" + + def group(self) -> str: + """Return the algorithm group name.""" + return "Contact Extractors" + + def groupId(self) -> str: + """Return the algorithm group ID.""" + return "Contact_Extractors" + + 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( + QgsProcessingParameterField( + 'UNIT_NAME_FIELD', + 'Unit Name Field', + parentLayerParameterName=self.INPUT_GEOLOGY, + type=QgsProcessingParameterField.String, + defaultValue='unitname' + ) + ) + + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_FAULTS, + "FAULTS", + [QgsProcessing.TypeVectorLine], + optional=True, + ) + ) + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_STRATI_COLUMN, + "Stratigraphic Order", + [QgsProcessing.TypeVector], + defaultValue='formation', + ) + ) + ignore_settings = QgsSettings() + last_ignore_units = ignore_settings.value("m2l/ignore_units", "") + self.addParameter( + QgsProcessingParameterMatrix( + self.INPUT_IGNORE_UNITS, + "Unit(s) to ignore", + headers=["Unit"], + defaultValue=last_ignore_units, + optional=True + ) + ) + + self.addParameter( + QgsProcessingParameterFeatureSink( + self.OUTPUT, + "Basal Contacts", + ) + ) + + self.addParameter( + QgsProcessingParameterFeatureSink( + "ALL_CONTACTS", + "All Contacts", + ) + ) + + def processAlgorithm( + self, + parameters: dict[str, Any], + context: QgsProcessingContext, + feedback: QgsProcessingFeedback, + ) -> dict[str, Any]: + + feedback.pushInfo("Loading data...") + geology = self.parameterAsVectorLayer(parameters, self.INPUT_GEOLOGY, context) + faults = self.parameterAsVectorLayer(parameters, self.INPUT_FAULTS, context) + strati_column = self.parameterAsSource(parameters, self.INPUT_STRATI_COLUMN, context) + ignore_units = self.parameterAsMatrix(parameters, self.INPUT_IGNORE_UNITS, context) + + + if isinstance(strati_column, QgsProcessingParameterMapLayer) : + raise QgsProcessingException("Invalid stratigraphic column layer") + + elif strati_column is not None: + # extract unit names from strati_column + field_name = "unit_name" + strati_order = [f[field_name] for f in strati_column.getFeatures()] + + if not ignore_units or all(isinstance(unit, str) and not unit.strip() for unit in ignore_units): + feedback.pushInfo("no units to ignore specified") + + ignore_settings = QgsSettings() + ignore_settings.setValue("m2l/ignore_units", ignore_units) + + unit_name_field = self.parameterAsString(parameters, 'UNIT_NAME_FIELD', context) + + geology = qgsLayerToGeoDataFrame(geology) + if unit_name_field and unit_name_field in geology.columns: + mask = ~geology[unit_name_field].astype(str).str.strip().isin(ignore_units) + geology = geology[mask].reset_index(drop=True) + feedback.pushInfo(f"filtered by unit name field: {unit_name_field}") + else: + feedback.pushInfo(f"no unit name field found: {unit_name_field}") + + faults = qgsLayerToGeoDataFrame(faults) if faults else None + if unit_name_field != 'UNITNAME' and unit_name_field in geology.columns: + geology = geology.rename(columns={unit_name_field: 'UNITNAME'}) + + feedback.pushInfo("Extracting Basal Contacts...") + contact_extractor = ContactExtractor(geology, faults) + all_contacts = contact_extractor.extract_all_contacts() + basal_contacts = contact_extractor.extract_basal_contacts(strati_order) + + feedback.pushInfo("Exporting Basal Contacts Layer...") + basal_contacts = GeoDataFrameToQgsLayer( + self, + basal_contacts, + parameters=parameters, + context=context, + output_key=self.OUTPUT, + feedback=feedback, + ) + contacts_layer = GeoDataFrameToQgsLayer( + self, + all_contacts, + parameters=parameters, + context=context, + output_key=self.ALL_CONTACTS, + feedback=feedback, + ) + return {self.OUTPUT: basal_contacts, self.ALL_CONTACTS: contacts_layer} + + def createInstance(self) -> QgsProcessingAlgorithm: + """Create a new instance of the algorithm.""" + return self.__class__() # BasalContactsAlgorithm() diff --git a/loopstructural/processing/algorithms/sampler.py b/loopstructural/processing/algorithms/sampler.py new file mode 100644 index 0000000..66191c3 --- /dev/null +++ b/loopstructural/processing/algorithms/sampler.py @@ -0,0 +1,225 @@ +""" +*************************************************************************** +* * +* 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, QVariant +from osgeo import gdal +import pandas as pd + +# QGIS imports +from qgis.core import ( + QgsFeatureSink, + QgsProcessing, + QgsProcessingAlgorithm, + QgsProcessingContext, + QgsProcessingException, + QgsProcessingFeedback, + QgsProcessingParameterFeatureSink, + QgsProcessingParameterFeatureSource, + QgsProcessingParameterRasterLayer, + QgsProcessingParameterEnum, + QgsProcessingParameterNumber, + QgsFields, + QgsField, + QgsFeature, + QgsGeometry, + QgsPointXY, + QgsVectorLayer, + QgsWkbTypes, + QgsCoordinateReferenceSystem +) +# Internal imports +from ...main.vectorLayerWrapper import qgsLayerToGeoDataFrame +from 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 "Spacing-Decimator Samplers" + + def group(self) -> str: + """Return the algorithm group name.""" + return "Samplers" + + def groupId(self) -> str: + """Return the algorithm group ID.""" + return "Samplers" + + def initAlgorithm(self, config: Optional[dict[str, Any]] = None) -> None: + """Initialize the algorithm parameters.""" + + + self.addParameter( + QgsProcessingParameterEnum( + self.INPUT_SAMPLER_TYPE, + "SAMPLER_TYPE", + ["Decimator (Point Geometry Data)", + "Spacing (Line Geometry Data)"], + defaultValue=0 + ) + ) + + self.addParameter( + QgsProcessingParameterRasterLayer( + self.INPUT_DTM, + "DTM", + [QgsProcessing.TypeRaster], + optional=True, + ) + ) + + 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 (Point Geometry Data)", + QgsProcessingParameterNumber.Integer, + defaultValue=1, + optional=True, + ) + ) + + self.addParameter( + QgsProcessingParameterNumber( + self.INPUT_SPACING, + "SPACING (Line Geometry Data)", + QgsProcessingParameterNumber.Double, + defaultValue=200.0, + optional=True, + ) + ) + + self.addParameter( + QgsProcessingParameterFeatureSink( + self.OUTPUT, + "Sampled Points", + ) + ) + + def processAlgorithm( + self, + parameters: dict[str, Any], + context: QgsProcessingContext, + feedback: QgsProcessingFeedback, + ) -> dict[str, Any]: + + dtm = self.parameterAsRasterLayer(parameters, self.INPUT_DTM, context) + geology = self.parameterAsVectorLayer(parameters, self.INPUT_GEOLOGY, context) + spatial_data = self.parameterAsVectorLayer(parameters, self.INPUT_SPATIAL_DATA, context) + decimation = self.parameterAsInt(parameters, self.INPUT_DECIMATION, context) + spacing = self.parameterAsDouble(parameters, self.INPUT_SPACING, context) + sampler_type_index = self.parameterAsEnum(parameters, self.INPUT_SAMPLER_TYPE, context) + sampler_type = ["Decimator", "Spacing"][sampler_type_index] + + if spatial_data is None: + raise QgsProcessingException("Spatial data is required") + + if sampler_type == "Decimator": + if geology is None: + raise QgsProcessingException("Geology is required") + if dtm is None: + raise QgsProcessingException("DTM is required") + + # Convert geology layers to GeoDataFrames + geology = qgsLayerToGeoDataFrame(geology) + spatial_data_gdf = qgsLayerToGeoDataFrame(spatial_data) + dtm_gdal = gdal.Open(dtm.source()) if dtm is not None and dtm.isValid() else None + + if sampler_type == "Decimator": + feedback.pushInfo("Sampling...") + sampler = SamplerDecimator(decimation=decimation, dtm_data=dtm_gdal, geology_data=geology) + samples = sampler.sample(spatial_data_gdf) + + if sampler_type == "Spacing": + feedback.pushInfo("Sampling...") + sampler = SamplerSpacing(spacing=spacing, dtm_data=dtm_gdal, geology_data=geology) + samples = sampler.sample(spatial_data_gdf) + + fields = QgsFields() + fields.append(QgsField("ID", QVariant.String)) + fields.append(QgsField("X", QVariant.Double)) + fields.append(QgsField("Y", QVariant.Double)) + fields.append(QgsField("Z", QVariant.Double)) + fields.append(QgsField("featureId", QVariant.String)) + + crs = None + if spatial_data_gdf is not None and spatial_data_gdf.crs is not None: + crs = QgsCoordinateReferenceSystem.fromWkt(spatial_data_gdf.crs.to_wkt()) + + sink, dest_id = self.parameterAsSink( + parameters, + self.OUTPUT, + context, + fields, + QgsWkbTypes.PointZ if 'Z' in (samples.columns if samples is not None else []) else QgsWkbTypes.Point, + crs + ) + + if samples is not None and not samples.empty: + for _index, row in samples.iterrows(): + feature = QgsFeature(fields) + + # decimator has z values + if 'Z' in samples.columns and pd.notna(row.get('Z')): + wkt = f"POINT Z ({row['X']} {row['Y']} {row['Z']})" + feature.setGeometry(QgsGeometry.fromWkt(wkt)) + else: + #spacing has no z values + feature.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(row['X'], row['Y']))) + + feature.setAttributes([ + str(row.get('ID', '')), + float(row.get('X', 0)), + float(row.get('Y', 0)), + float(row.get('Z', 0)) if pd.notna(row.get('Z')) else 0.0, + str(row.get('featureId', '')) + ]) + + sink.addFeature(feature) + + 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/loopstructural/processing/algorithms/sorter.py b/loopstructural/processing/algorithms/sorter.py new file mode 100644 index 0000000..55a179c --- /dev/null +++ b/loopstructural/processing/algorithms/sorter.py @@ -0,0 +1,418 @@ +from typing import Any, Optional +from osgeo import gdal +import pandas as pd +import json + +from PyQt5.QtCore import QVariant +from qgis import processing +from qgis.core import ( + QgsFeatureSink, + QgsFields, + QgsField, + QgsFeature, + QgsGeometry, + QgsRasterLayer, + QgsProcessing, + QgsProcessingAlgorithm, + QgsProcessingContext, + QgsProcessingException, + QgsProcessingFeedback, + QgsProcessingParameterEnum, + QgsProcessingParameterFileDestination, + QgsProcessingParameterFeatureSink, + QgsProcessingParameterFeatureSource, + QgsProcessingParameterField, + QgsProcessingParameterRasterLayer, + QgsProcessingParameterMatrix, + QgsCoordinateReferenceSystem, + QgsVectorLayer, + QgsWkbTypes, + QgsSettings +) + +# ──────────────────────────────────────────────── +# map2loop sorters +# ──────────────────────────────────────────────── +from map2loop.sorter import ( + SorterAlpha, + SorterAgeBased, + SorterMaximiseContacts, + SorterObservationProjections, + SorterUseNetworkX, + SorterUseHint, # kept for backwards compatibility +) +from ...main.vectorLayerWrapper import qgsLayerToGeoDataFrame, qvariantToFloat + +# 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. + """ + METHOD = "METHOD" + INPUT_GEOLOGY = "INPUT_GEOLOGY" + INPUT_STRUCTURE = "INPUT_STRUCTURE" + INPUT_DTM = "INPUT_DTM" + INPUT_STRATI_COLUMN = "INPUT_STRATI_COLUMN" + SORTING_ALGORITHM = "SORTING_ALGORITHM" + OUTPUT = "OUTPUT" + CONTACTS_LAYER = "CONTACTS_LAYER" + + # ---------------------------------------------------------- + # Metadata + # ---------------------------------------------------------- + def name(self) -> str: + return "loop_sorter" + + def displayName(self) -> str: + return "Automatic Stratigraphic Column" + + def group(self) -> str: + return "Stratigraphy" + + def groupId(self) -> str: + return "Stratigraphy_Column" + + def updateParameters(self, parameters): + selected_method = parameters.get(self.METHOD, 0) + selected_algorithm = parameters.get(self.SORTING_ALGORITHM, 0) + + if selected_method == 0: # User-Defined selected + self.parameterDefinition(self.INPUT_STRATI_COLUMN).setMetadata({'widget_wrapper': {'visible': True}}) + self.parameterDefinition(self.SORTING_ALGORITHM).setMetadata({'widget_wrapper': {'visible': False}}) + self.parameterDefinition(self.INPUT_GEOLOGY).setMetadata({'widget_wrapper': {'visible': False}}) + else: # Automatic selected + self.parameterDefinition(self.INPUT_GEOLOGY).setMetadata({'widget_wrapper': {'visible': True}}) + self.parameterDefinition(self.SORTING_ALGORITHM).setMetadata({'widget_wrapper': {'visible': True}}) + self.parameterDefinition(self.INPUT_STRATI_COLUMN).setMetadata({'widget_wrapper': {'visible': False}}) + + # observation projects + is_observation_projections = selected_algorithm == 5 + self.parameterDefinition(self.INPUT_STRUCTURE).setMetadata({'widget_wrapper': {'visible': is_observation_projections}}) + self.parameterDefinition(self.INPUT_DTM).setMetadata({'widget_wrapper': {'visible': is_observation_projections}}) + self.parameterDefinition('DIP_FIELD').setMetadata({'widget_wrapper': {'visible': is_observation_projections}}) + self.parameterDefinition('DIPDIR_FIELD').setMetadata({'widget_wrapper': {'visible': is_observation_projections}}) + self.parameterDefinition('ORIENTATION_TYPE').setMetadata({'widget_wrapper': {'visible': is_observation_projections}}) + + return super().updateParameters(parameters) + + # ---------------------------------------------------------- + # Parameters + # ---------------------------------------------------------- + def initAlgorithm(self, config: Optional[dict[str, Any]] = None) -> None: + + # enum so the user can pick the strategy from a dropdown + self.addParameter( + QgsProcessingParameterEnum( + self.SORTING_ALGORITHM, + "Sorting strategy", + options=list(SORTER_LIST.keys()), + defaultValue="Observation projections", # Age-based is safest default + ) + ) + + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_GEOLOGY, + "Geology polygons", + [QgsProcessing.TypeVectorPolygon], + optional=True + ) + ) + + self.addParameter( + QgsProcessingParameterField( + 'UNIT_NAME_FIELD', + 'Unit Name Field', + parentLayerParameterName=self.INPUT_GEOLOGY, + type=QgsProcessingParameterField.Any, + defaultValue='UNITNAME', + optional=True + ) + ) + + self.addParameter( + QgsProcessingParameterField( + 'MIN_AGE_FIELD', + 'Minimum Age Field', + parentLayerParameterName=self.INPUT_GEOLOGY, + type=QgsProcessingParameterField.Any, + defaultValue='MIN_AGE', + optional=True + ) + ) + + self.addParameter( + QgsProcessingParameterField( + 'MAX_AGE_FIELD', + 'Maximum Age Field', + parentLayerParameterName=self.INPUT_GEOLOGY, + type=QgsProcessingParameterField.Any, + defaultValue='MAX_AGE', + optional=True + ) + ) + + self.addParameter( + QgsProcessingParameterField( + 'GROUP_FIELD', + 'Group Field', + parentLayerParameterName=self.INPUT_GEOLOGY, + type=QgsProcessingParameterField.Any, + defaultValue='GROUP', + optional=True + ) + ) + + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_STRUCTURE, + "Structure", + [QgsProcessing.TypeVectorPoint], + optional=True, + ) + ) + + self.addParameter( + QgsProcessingParameterField( + 'DIP_FIELD', + 'Dip Field', + parentLayerParameterName=self.INPUT_STRUCTURE, + type=QgsProcessingParameterField.Any, + defaultValue='DIP', + optional=True + ) + ) + self.addParameter( + QgsProcessingParameterField( + 'DIPDIR_FIELD', + 'Dip Direction Field', + parentLayerParameterName=self.INPUT_STRUCTURE, + type=QgsProcessingParameterField.Any, + defaultValue='DIPDIR', + optional=True + ) + ) + + self.addParameter( + QgsProcessingParameterEnum( + 'ORIENTATION_TYPE', + 'Orientation Type', + options=['','Dip Direction', 'Strike'], + defaultValue=0 + ) + ) + + self.addParameter( + QgsProcessingParameterRasterLayer( + self.INPUT_DTM, + "DTM", + optional=True, + ) + ) + + self.addParameter( + QgsProcessingParameterFeatureSource( + "CONTACTS_LAYER", + "Contacts Layer", + [QgsProcessing.TypeVectorLine], + optional=False, + ) + ) + + self.addParameter( + QgsProcessingParameterFeatureSink( + self.OUTPUT, + "Stratigraphic column", + ) + ) + + self.addParameter( + QgsProcessingParameterFileDestination( + "JSON_OUTPUT", + "Stratigraphic column json", + fileFilter="JSON files (*.json)" + ) + ) + + # ---------------------------------------------------------- + # Core + # ---------------------------------------------------------- + def processAlgorithm( + self, + parameters: dict[str, Any], + context: QgsProcessingContext, + feedback: QgsProcessingFeedback, + ) -> dict[str, Any]: + + algo_index: int = self.parameterAsEnum(parameters, self.SORTING_ALGORITHM, context) + sorter_cls = list(SORTER_LIST.values())[algo_index] + contacts_layer = self.parameterAsVectorLayer(parameters, self.CONTACTS_LAYER, context) + in_layer = self.parameterAsVectorLayer(parameters, self.INPUT_GEOLOGY, context) + output_file = self.parameterAsFileOutput(parameters, 'JSON_OUTPUT', context) + + units_df, relationships_df, contacts_df= build_input_frames(in_layer,contacts_layer, feedback,parameters) + + if sorter_cls == SorterObservationProjections: + geology_gdf = qgsLayerToGeoDataFrame(in_layer) + structure = self.parameterAsVectorLayer(parameters, self.INPUT_STRUCTURE, context) + dtm = self.parameterAsRasterLayer(parameters, self.INPUT_DTM, context) + if geology_gdf is None or geology_gdf.empty or not structure or not structure.isValid() or not dtm or not dtm.isValid(): + raise QgsProcessingException("Structure and DTM layer are required for observation projections") + + structure_gdf = qgsLayerToGeoDataFrame(structure) if structure else None + dtm_gdal = gdal.Open(dtm.source()) if dtm is not None and dtm.isValid() else None + + unit_name_field = parameters.get('UNIT_NAME_FIELD', 'UNITNAME') if parameters else 'UNITNAME' + if unit_name_field != 'UNITNAME' and unit_name_field in geology_gdf.columns: + geology_gdf = geology_gdf.rename(columns={unit_name_field: 'UNITNAME'}) + + dip_field = parameters.get('DIP_FIELD', 'DIP') if parameters else 'DIP' + if not dip_field: + raise QgsProcessingException("Dip Field is required") + if dip_field != 'DIP' and dip_field in structure_gdf.columns: + structure_gdf = structure_gdf.rename(columns={dip_field: 'DIP'}) + orientation_type = self.parameterAsEnum(parameters, 'ORIENTATION_TYPE', context) + orientation_type_name = ['','Dip Direction', 'Strike'][orientation_type] + if not orientation_type_name: + raise QgsProcessingException("Orientation Type is required") + dipdir_field = parameters.get('DIPDIR_FIELD', 'DIPDIR') if parameters else 'DIPDIR' + if not dipdir_field: + raise QgsProcessingException("Dip Direction Field is required") + if dipdir_field in structure_gdf.columns: + if orientation_type_name == 'Strike': + structure_gdf['DIPDIR'] = structure_gdf[dipdir_field].apply( + lambda val: (val + 90.0) % 360.0 if pd.notnull(val) else val + ) + elif orientation_type_name == 'Dip Direction': + structure_gdf = structure_gdf.rename(columns={dipdir_field: 'DIPDIR'}) + order = sorter_cls().sort( + units_df, + relationships_df, + contacts_df, + geology_gdf, + structure_gdf, + dtm_gdal + ) + else: + order = sorter_cls().sort( + units_df, + relationships_df, + contacts_df + ) + + # 4 β–Ί write an in-memory table with the result + sink_fields = QgsFields() + sink_fields.append(QgsField("order", QVariant.Int)) + sink_fields.append(QgsField("unit_name", QVariant.String)) + + (sink, dest_id) = self.parameterAsSink( + parameters, + self.OUTPUT, + context, + sink_fields, + QgsWkbTypes.NoGeometry, + in_layer.sourceCrs() if in_layer else None, + ) + + for pos, name in enumerate(order, start=1): + f = QgsFeature(sink_fields) + f.setAttributes([pos, name]) + sink.addFeature(f, QgsFeatureSink.FastInsert) + try: + with open(output_file, 'w') as f: + json.dump(order, f) + except Exception as e: + with open(output_file, 'w') as f: + json.dump([], f) + + return {self.OUTPUT: dest_id, 'JSON_OUTPUT': output_file} + + # ---------------------------------------------------------- + def createInstance(self) -> QgsProcessingAlgorithm: + return StratigraphySorterAlgorithm() + +# ------------------------------------------------------------------------- +# Helper stub – you must replace with *your* conversion logic +# ------------------------------------------------------------------------- +def build_input_frames(layer: QgsVectorLayer,contacts_layer: QgsVectorLayer, feedback, parameters, user_defined_units=None) -> 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) + """ + + if user_defined_units: + units_record = [] + for i, row in enumerate(user_defined_units): + units_record.append( + dict( + layerId=i, + name=row[1], + minAge=row[2], + maxAge=row[3], + group=row[4] + ) + ) + units_df = pd.DataFrame.from_records(units_record) + else: + unit_name_field = parameters.get('UNIT_NAME_FIELD', 'UNITNAME') if parameters else 'UNITNAME' + min_age_field = parameters.get('MIN_AGE_FIELD', 'MIN_AGE') if parameters else 'MIN_AGE' + max_age_field = parameters.get('MAX_AGE_FIELD', 'MAX_AGE') if parameters else 'MAX_AGE' + group_field = parameters.get('GROUP_FIELD', 'GROUP') if parameters else 'GROUP' + + if not layer or not layer.isValid(): + raise QgsProcessingException("No geology layer provided") + if not unit_name_field: + raise QgsProcessingException("Unit Name Field is required") + if not min_age_field: + raise QgsProcessingException("Minimum Age Field is required") + if not max_age_field: + raise QgsProcessingException("Maximum Age Field is required") + if not group_field: + raise QgsProcessingException("Group Field is required") + + units_records = [] + for f in layer.getFeatures(): + units_records.append( + dict( + layerId=f.id(), + name=f[unit_name_field], + minAge=qvariantToFloat(f, min_age_field), + maxAge=qvariantToFloat(f, max_age_field), + group=f[group_field], + ) + ) + units_df = pd.DataFrame.from_records(units_records) + + feedback.pushInfo(f"Units β†’ {len(units_df)} records") + # map_data can be mocked if you only use Age-based sorter + + if not contacts_layer or not contacts_layer.isValid(): + raise QgsProcessingException("No contacts layer provided") + + contacts_df = qgsLayerToGeoDataFrame(contacts_layer) if contacts_layer else pd.DataFrame() + if not contacts_df.empty: + relationships_df = contacts_df.copy() + if 'length' in contacts_df.columns: + relationships_df = relationships_df.drop(columns=['length']) + if 'geometry' in contacts_df.columns: + relationships_df = relationships_df.drop(columns=['geometry']) + feedback.pushInfo(f"Contacts β†’ {len(contacts_df)} records") + feedback.pushInfo(f"Relationships β†’ {len(relationships_df)} records") + else: + relationships_df = pd.DataFrame() + + return units_df, relationships_df, contacts_df diff --git a/loopstructural/processing/algorithms/thickness_calculator.py b/loopstructural/processing/algorithms/thickness_calculator.py new file mode 100644 index 0000000..8d67dc5 --- /dev/null +++ b/loopstructural/processing/algorithms/thickness_calculator.py @@ -0,0 +1,377 @@ +""" +*************************************************************************** +* * +* 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 +import pandas as pd + +# QGIS imports +from qgis import processing +from qgis.core import ( + QgsFeatureSink, + QgsProcessing, + QgsProcessingAlgorithm, + QgsProcessingContext, + QgsProcessingException, + QgsProcessingFeedback, + QgsProcessingParameterFeatureSink, + QgsProcessingParameterFeatureSource, + QgsProcessingParameterEnum, + QgsProcessingParameterNumber, + QgsProcessingParameterField, + QgsProcessingParameterMatrix, + QgsSettings, + QgsProcessingParameterRasterLayer, +) +# Internal imports +from ...main.vectorLayerWrapper import ( + qgsLayerToGeoDataFrame, + GeoDataFrameToQgsLayer, + qgsLayerToDataFrame, + dataframeToQgsLayer, + qgsRasterToGdalDataset, + matrixToDict, + dataframeToQgsTable + ) +from 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_TYPE = 'BOUNDING_BOX_TYPE' + INPUT_BOUNDING_BOX = 'BOUNDING_BOX' + INPUT_MAX_LINE_LENGTH = 'MAX_LINE_LENGTH' + INPUT_STRATI_COLUMN = 'STRATIGRAPHIC_COLUMN' + INPUT_BASAL_CONTACTS = 'BASAL_CONTACTS' + INPUT_STRUCTURE_DATA = 'STRUCTURE_DATA' + INPUT_DIPDIR_FIELD = 'DIPDIR_FIELD' + INPUT_DIP_FIELD = 'DIP_FIELD' + INPUT_GEOLOGY = 'GEOLOGY' + INPUT_ORIENTATION_TYPE = 'ORIENTATION_TYPE' + INPUT_UNIT_NAME_FIELD = 'UNIT_NAME_FIELD' + INPUT_SAMPLED_CONTACTS = 'SAMPLED_CONTACTS' + INPUT_STRATIGRAPHIC_COLUMN_LAYER = 'STRATIGRAPHIC_COLUMN_LAYER' + + OUTPUT = "THICKNESS" + + def name(self) -> str: + """Return the algorithm name.""" + return "thickness_calculator" + + def displayName(self) -> str: + """Return the algorithm display name.""" + return "Thickness Calculator" + + def group(self) -> str: + """Return the algorithm group name.""" + return "Thickness Calculators" + + def groupId(self) -> str: + """Return the algorithm group ID.""" + return "Thickness_Calculators" + + 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, + defaultValue='InterpolatedStructure' + ) + ) + self.addParameter( + QgsProcessingParameterRasterLayer( + self.INPUT_DTM, + "DTM (InterpolatedStructure)", + [QgsProcessing.TypeRaster], + optional=True, + ) + ) + + self.addParameter( + QgsProcessingParameterEnum( + self.INPUT_BOUNDING_BOX_TYPE, + "Bounding Box Type", + options=['Extract from geology layer', 'User defined'], + allowMultiple=False, + defaultValue=1 + ) + ) + + bbox_settings = QgsSettings() + last_bbox = bbox_settings.value("m2l/bounding_box", "") + self.addParameter( + QgsProcessingParameterMatrix( + self.INPUT_BOUNDING_BOX, + description="Static Bounding Box", + headers=['minx','miny','maxx','maxy'], + numberRows=1, + defaultValue=last_bbox, + optional=True + ) + ) + + self.addParameter( + QgsProcessingParameterNumber( + self.INPUT_MAX_LINE_LENGTH, + "Max Line Length", + minValue=0, + defaultValue=1000 + ) + ) + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_BASAL_CONTACTS, + "Basal Contacts", + [QgsProcessing.TypeVectorLine], + defaultValue='Basal Contacts', + ) + ) + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_GEOLOGY, + "GEOLOGY", + [QgsProcessing.TypeVectorPolygon], + ) + ) + + self.addParameter( + QgsProcessingParameterField( + 'UNIT_NAME_FIELD', + 'Unit Name Field e.g. Formation', + parentLayerParameterName=self.INPUT_GEOLOGY, + type=QgsProcessingParameterField.String, + defaultValue='Formation' + ) + ) + + self.addParameter( + QgsProcessingParameterFeatureSource( + 'STRATIGRAPHIC_COLUMN_LAYER', + 'Stratigraphic Column Layer (from sorter)', + [QgsProcessing.TypeVector], + optional=True + ) + ) + + strati_settings = QgsSettings() + last_strati_column = strati_settings.value("m2l/strati_column", "") + self.addParameter( + QgsProcessingParameterMatrix( + name=self.INPUT_STRATI_COLUMN, + description="Stratigraphic Order", + headers=["Unit"], + numberRows=0, + defaultValue=last_strati_column, + optional=True + ) + ) + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_SAMPLED_CONTACTS, + "Sampled Contacts", + [QgsProcessing.TypeVectorPoint], + ) + ) + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_STRUCTURE_DATA, + "Orientation Data", + [QgsProcessing.TypeVectorPoint], + ) + ) + self.addParameter( + QgsProcessingParameterEnum( + self.INPUT_ORIENTATION_TYPE, + 'Orientation Type', + options=['Dip Direction', 'Strike'], + defaultValue=0 # Default to Dip Direction + ) + ) + self.addParameter( + QgsProcessingParameterField( + self.INPUT_DIPDIR_FIELD, + "Dip Direction Column", + parentLayerParameterName=self.INPUT_STRUCTURE_DATA, + type=QgsProcessingParameterField.Numeric, + defaultValue='DIPDIR' + ) + ) + self.addParameter( + QgsProcessingParameterField( + self.INPUT_DIP_FIELD, + "Dip Column", + parentLayerParameterName=self.INPUT_STRUCTURE_DATA, + type=QgsProcessingParameterField.Numeric, + defaultValue='DIP' + ) + ) + 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_index = self.parameterAsEnum(parameters, self.INPUT_THICKNESS_CALCULATOR_TYPE, context) + thickness_type = ['InterpolatedStructure', 'StructuralPoint'][thickness_type_index] + dtm_data = self.parameterAsRasterLayer(parameters, self.INPUT_DTM, context) + bounding_box_type = self.parameterAsEnum(parameters, self.INPUT_BOUNDING_BOX_TYPE, context) + max_line_length = self.parameterAsSource(parameters, self.INPUT_MAX_LINE_LENGTH, context) + basal_contacts = self.parameterAsSource(parameters, self.INPUT_BASAL_CONTACTS, context) + geology_data = self.parameterAsSource(parameters, self.INPUT_GEOLOGY, context) + structure_data = self.parameterAsSource(parameters, self.INPUT_STRUCTURE_DATA, context) + orientation_type = self.parameterAsEnum(parameters, self.INPUT_ORIENTATION_TYPE, context) + is_strike = (orientation_type == 1) + structure_dipdir_field = self.parameterAsString(parameters, self.INPUT_DIPDIR_FIELD, context) + structure_dip_field = self.parameterAsString(parameters, self.INPUT_DIP_FIELD, context) + sampled_contacts = self.parameterAsSource(parameters, self.INPUT_SAMPLED_CONTACTS, context) + unit_name_field = self.parameterAsString(parameters, self.INPUT_UNIT_NAME_FIELD, context) + + if bounding_box_type == 0: + geology_layer = self.parameterAsVectorLayer(parameters, self.INPUT_GEOLOGY, context) + extent = geology_layer.extent() + bounding_box = { + 'minx': extent.xMinimum(), + 'miny': extent.yMinimum(), + 'maxx': extent.xMaximum(), + 'maxy': extent.yMaximum() + } + feedback.pushInfo("Using bounding box from geology layer") + else: + static_bbox_matrix = self.parameterAsMatrix(parameters, self.INPUT_BOUNDING_BOX, context) + if not static_bbox_matrix or len(static_bbox_matrix) == 0: + raise QgsProcessingException("Bounding box is required") + + bounding_box = matrixToDict(static_bbox_matrix) + + bbox_settings = QgsSettings() + bbox_settings.setValue("m2l/bounding_box", static_bbox_matrix) + feedback.pushInfo("Using bounding box from user input") + + stratigraphic_column_source = self.parameterAsSource(parameters, self.INPUT_STRATIGRAPHIC_COLUMN_LAYER, context) + stratigraphic_order = [] + if stratigraphic_column_source is not None: + ordered_pairs =[] + for feature in stratigraphic_column_source.getFeatures(): + order = feature.attribute('order') + unit_name = feature.attribute('unit_name') + ordered_pairs.append((order, unit_name)) + ordered_pairs.sort(key=lambda x: x[0]) + stratigraphic_order = [pair[1] for pair in ordered_pairs] + feedback.pushInfo(f"DEBUG: parameterAsVectorLayer Stratigraphic order: {stratigraphic_order}") + else: + matrix_stratigraphic_order = self.parameterAsMatrix(parameters, self.INPUT_STRATI_COLUMN, context) + if matrix_stratigraphic_order: + stratigraphic_order = [str(row) for row in matrix_stratigraphic_order if row] + else: + raise QgsProcessingException("Stratigraphic column layer is required") + if stratigraphic_order: + matrix = [[unit] for unit in stratigraphic_order] + strati_column_settings = QgsSettings() + strati_column_settings.setValue('m2l/strati_column', matrix) + # convert layers to dataframe or geodataframe + units = qgsLayerToDataFrame(geology_data) + geology_data = qgsLayerToGeoDataFrame(geology_data) + basal_contacts = qgsLayerToGeoDataFrame(basal_contacts) + structure_data = qgsLayerToDataFrame(structure_data) + rename_map = {} + missing_fields = [] + if unit_name_field != 'UNITNAME' and unit_name_field in geology_data.columns: + geology_data = geology_data.rename(columns={unit_name_field: 'UNITNAME'}) + units_unique = units.drop_duplicates(subset=[unit_name_field]).reset_index(drop=True) + units = pd.DataFrame({'name': units_unique[unit_name_field]}) + if structure_data is not None: + if structure_dipdir_field: + if structure_dipdir_field in structure_data.columns: + rename_map[structure_dipdir_field] = 'DIPDIR' + else: + missing_fields.append(structure_dipdir_field) + if structure_dip_field: + if structure_dip_field in structure_data.columns: + rename_map[structure_dip_field] = 'DIP' + else: + missing_fields.append(structure_dip_field) + if missing_fields: + raise QgsProcessingException( + f"Orientation data missing required field(s): {', '.join(missing_fields)}" + ) + if rename_map: + structure_data = structure_data.rename(columns=rename_map) + + sampled_contacts = qgsLayerToDataFrame(sampled_contacts) + sampled_contacts['X'] = sampled_contacts['X'].astype(float) + sampled_contacts['Y'] = sampled_contacts['Y'].astype(float) + sampled_contacts['Z'] = sampled_contacts['Z'].astype(float) + dtm_data = qgsRasterToGdalDataset(dtm_data) + if thickness_type == "InterpolatedStructure": + thickness_calculator = InterpolatedStructure( + dtm_data=dtm_data, + bounding_box=bounding_box, + is_strike=is_strike + ) + thicknesses = 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, + is_strike=is_strike + ) + thicknesses =thickness_calculator.compute( + units, + stratigraphic_order, + basal_contacts, + structure_data, + geology_data, + sampled_contacts + ) + + thicknesses = thicknesses[ + ["name","ThicknessMean","ThicknessMedian", "ThicknessStdDev"] + ].copy() + + feedback.pushInfo("Exporting Thickness Table...") + thicknesses = dataframeToQgsTable( + self, + thicknesses, + parameters=parameters, + context=context, + feedback=feedback, + param_name=self.OUTPUT + ) + + return {self.OUTPUT: thicknesses[1]} + + def createInstance(self) -> QgsProcessingAlgorithm: + """Create a new instance of the algorithm.""" + return self.__class__() # ThicknessCalculatorAlgorithm() diff --git a/loopstructural/processing/algorithms/user_defined_sorter.py b/loopstructural/processing/algorithms/user_defined_sorter.py new file mode 100644 index 0000000..23857a3 --- /dev/null +++ b/loopstructural/processing/algorithms/user_defined_sorter.py @@ -0,0 +1,112 @@ +from typing import Any, Optional +from osgeo import gdal +import numpy as np +import json + +from PyQt5.QtCore import QVariant +from qgis import processing +from qgis.core import ( + QgsFeatureSink, + QgsFields, + QgsField, + QgsFeature, + QgsGeometry, + QgsRasterLayer, + QgsProcessing, + QgsProcessingAlgorithm, + QgsProcessingContext, + QgsProcessingException, + QgsProcessingFeedback, + QgsProcessingParameterEnum, + QgsProcessingParameterFileDestination, + QgsProcessingParameterFeatureSink, + QgsProcessingParameterFeatureSource, + QgsProcessingParameterField, + QgsProcessingParameterRasterLayer, + QgsProcessingParameterMatrix, + QgsCoordinateReferenceSystem, + QgsVectorLayer, + QgsWkbTypes, + QgsSettings +) + + +from qgis.core import ( + QgsFields, QgsField, QgsFeature, QgsFeatureSink, QgsWkbTypes, + QgsCoordinateReferenceSystem, QgsProcessingAlgorithm, QgsProcessingContext, + QgsProcessingFeedback, QgsProcessingParameterFeatureSink, QgsProcessingParameterMatrix, + QgsSettings +) +from PyQt5.QtCore import QVariant +import numpy as np + +class UserDefinedStratigraphyAlgorithm(QgsProcessingAlgorithm): + INPUT_STRATI_COLUMN = "INPUT_STRATI_COLUMN" + OUTPUT = "OUTPUT" + + def name(self): return "loop_sorter_2" + def displayName(self): return "User-Defined Stratigraphic Column" + def group(self): return "Stratigraphy" + def groupId(self): return "Stratigraphy_Column" + + def initAlgorithm(self, config=None): + strati_settings = QgsSettings() + last_strati_column = strati_settings.value("m2l/strati_column", "") + self.addParameter( + QgsProcessingParameterMatrix( + name=self.INPUT_STRATI_COLUMN, + description="Stratigraphic Order", + headers=["Unit"], + numberRows=0, + defaultValue=last_strati_column + ) + ) + self.addParameter( + QgsProcessingParameterFeatureSink( + self.OUTPUT, + "Stratigraphic column", + ) + ) + + def processAlgorithm(self, parameters, context, feedback): + # 1) Read the matrix; it may be a list of lists (rows) or a flat list depending on input source. + matrix = self.parameterAsMatrix(parameters, self.INPUT_STRATI_COLUMN, context) + + # Normalize to a list of unit strings (one column: "Unit") + units = [] + for row in matrix: + if isinstance(row, (list, tuple)): + unit = row[0] if row else "" + else: + unit = row + unit = (unit or "").strip() + if unit: # skip empty rows to avoid writing "" into fields + units.append(unit) + + # 2) Build sequential order (1-based), cast to native int + order_vals = [int(i) for i in (np.arange(len(units)) + 1)] + + # 3) Prepare sink + sink_fields = QgsFields() + sink_fields.append(QgsField("order", QVariant.Int)) # or QVariant.LongLong + sink_fields.append(QgsField("unit_name", QVariant.String)) + + crs = context.project().crs() if context and context.project() else QgsCoordinateReferenceSystem() + sink, dest_id = self.parameterAsSink( + parameters, self.OUTPUT, context, + sink_fields, QgsWkbTypes.NoGeometry, crs + ) + + # 4) Insert features + for pos, unit_name in zip(order_vals, units): + f = QgsFeature(sink_fields) + # Ensure correct types: int for "order", str for "unit_name" + f.setAttributes([int(pos), str(unit_name)]) + ok = sink.addFeature(f, QgsFeatureSink.FastInsert) + if not ok: + feedback.reportError(f"Failed to add feature for unit '{unit_name}' (order={pos}).") + + return {self.OUTPUT: dest_id} + + def createInstance(self): + return __class__() diff --git a/loopstructural/processing/provider.py b/loopstructural/processing/provider.py index 0b970cb..2e1fef0 100644 --- a/loopstructural/processing/provider.py +++ b/loopstructural/processing/provider.py @@ -8,7 +8,19 @@ from qgis.PyQt.QtGui import QIcon # project -from loopstructural.__about__ import __icon_path__, __title__, __version__ +from loopstructural.__about__ import ( + __icon_path__, + __title__, + __version__, +) + +from .algorithms import ( + BasalContactsAlgorithm, + SamplerAlgorithm, + StratigraphySorterAlgorithm, + ThicknessCalculatorAlgorithm, + UserDefinedStratigraphyAlgorithm, +) # ############################################################################ # ########## Classes ############### @@ -40,51 +52,77 @@ def name(self) -> str: ------- str Short, localised provider name. + + """ + return __title__ + + +class Map2LoopProvider(QgsProcessingProvider): + """Processing provider class.""" + + def loadAlgorithms(self): + """Loads all algorithms belonging to this provider.""" + self.addAlgorithm(BasalContactsAlgorithm()) + self.addAlgorithm(StratigraphySorterAlgorithm()) + self.addAlgorithm(UserDefinedStratigraphyAlgorithm()) + self.addAlgorithm(ThicknessCalculatorAlgorithm()) + self.addAlgorithm(SamplerAlgorithm()) + + def id(self) -> str: + """Unique provider id, used for identifying it. This string should be unique, \ + short, character only string, eg "qgis" or "gdal". \ + This string should not be localised. + + :return: provider ID + :rtype: str + """ + return "plugin_map2loop" + + def name(self) -> str: + """Returns the provider name, which is used to describe the provider + within the GUI. This string should be short (e.g. "Lastools") and localised. + + :return: provider name + :rtype: str """ return __title__ def longName(self) -> str: - """Return longer provider name (may include version information). + """Longer version of the provider name, which can include + extra details such as version numbers. E.g. "Lastools LIDAR tools". This string should be localised. The default + implementation returns the same string as name(). - Returns - ------- - str - Localised long name for display in the GUI. + :return: provider long name + :rtype: str """ return self.tr("{} - Tools".format(__title__)) def icon(self) -> QIcon: - """Return icon used for the provider inside the Processing toolbox menu. + """QIcon used for your provider inside the Processing toolbox menu. - Returns - ------- - QIcon - Icon for the provider. + :return: provider icon + :rtype: QIcon """ return QIcon(str(__icon_path__)) def tr(self, message: str) -> str: - """Translate a string using Qt translation API. + """Get the translation for a string using Qt translation API. - Parameters - ---------- - message : str - String to be translated. + :param message: String for translation. + :type message: str, QString - Returns - ------- - str - Translated version of message. + :returns: Translated version of message. + :rtype: str """ # noinspection PyTypeChecker,PyArgumentList,PyCallByClass return QCoreApplication.translate(self.__class__.__name__, message) def versionInfo(self) -> str: - """Return provider version information. + """Version information for the provider, or an empty string if this is not \ + applicable (e.g. for inbuilt Processing providers). For plugin based providers, \ + this should return the plugin’s version identifier. - Returns - ------- - str - Version string for the provider (plugin version). + :return: version + :rtype: str """ return __version__ diff --git a/loopstructural/resources/i18n/plugin_map2loop_en.ts b/loopstructural/resources/i18n/plugin_map2loop_en.ts new file mode 100644 index 0000000..4c23f0c --- /dev/null +++ b/loopstructural/resources/i18n/plugin_map2loop_en.ts @@ -0,0 +1,4 @@ + + + + diff --git a/loopstructural/toolbelt/env_var_parser.py b/loopstructural/toolbelt/env_var_parser.py new file mode 100644 index 0000000..2baa665 --- /dev/null +++ b/loopstructural/toolbelt/env_var_parser.py @@ -0,0 +1,60 @@ +import os +from typing import Type, TypeVar + +T = TypeVar("T") + + +class EnvVarParser: + """Utility class to retrieve and convert environment variables.""" + + @staticmethod + def get_env_var(name: str, default: T) -> T: + """Retrieves an environment variable and converts it based on the default value type. + + Args: + name (str): The environment variable name. + default (T): The default value, used to infer the expected type. + + Returns: + T: The converted value, matching the type of `default`. + """ + value = os.getenv(name) + if value is None: + return ( + default # Return the default value if the environment variable is not + ) + + # Otherwise, treat it as a single value + return EnvVarParser._convert_single(value, type(default), default) + + @staticmethod + def _convert_single(value: str, expected_type: Type[T], default: T) -> T: + """Converts a string into a single value of the expected type.""" + try: + if expected_type is int: + return int(value) + elif expected_type is float: + return float(value) + elif expected_type is bool: + return EnvVarParser._convert_bool(value, default) + elif expected_type is str: + return value # String value + except ValueError: + return default # Return default value in case of conversion failure + + raise TypeError( + f"Unsupported type: {expected_type}. Value definition from environment variable is not possible." + ) + + @staticmethod + def _convert_bool(value: str, default: bool) -> bool: + """Converts a string into a boolean, handling explicit True/False values.""" + true_values = {"1", "true", "yes", "on"} + false_values = {"0", "false", "no", "off"} + value_lower = value.lower() + + if value_lower in true_values: + return True + elif value_lower in false_values: + return False + return default # Return default value if conversion fails diff --git a/loopstructural/toolbelt/log_handler.py b/loopstructural/toolbelt/log_handler.py index 31001a0..0e97360 100644 --- a/loopstructural/toolbelt/log_handler.py +++ b/loopstructural/toolbelt/log_handler.py @@ -9,10 +9,10 @@ # standard library import logging from functools import partial -from typing import Callable +from typing import Callable, Literal, Optional, Union # PyQGIS -from qgis.core import QgsMessageLog, QgsMessageOutput +from qgis.core import Qgis, QgsMessageLog, QgsMessageOutput from qgis.gui import QgsMessageBar from qgis.PyQt.QtWidgets import QPushButton, QWidget from qgis.utils import iface diff --git a/requirements/development.txt b/requirements/development.txt index 7942959..703f00a 100644 --- a/requirements/development.txt +++ b/requirements/development.txt @@ -5,4 +5,4 @@ black isort>=5.13 -pre-commit>=3,<5 +pre-commit>=4,<5 diff --git a/requirements/testing.txt b/requirements/testing.txt index 53938b7..954e649 100644 --- a/requirements/testing.txt +++ b/requirements/testing.txt @@ -7,3 +7,5 @@ pyvista>=0.4 #3D visualisation/vtk meshing LoopStructural>=1.6.6 #3D geological modelling geoh5py>=0.9 #saving models to hdf5 dill>=0.3.9 #saving models to python +shapely +geopandas diff --git a/scripts/__init__.py b/scripts/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/scripts/generate_translation_profile.py b/scripts/generate_translation_profile.py new file mode 100644 index 0000000..5536d3e --- /dev/null +++ b/scripts/generate_translation_profile.py @@ -0,0 +1,86 @@ +#! python3 + +"""Script to generate the translation profile for a PyQt project, listing Python, +forms (ui) and targetted translation files. This script is intended to be run from the +root of the project, and will create a file named `plugin_translation.pro` in the +`oslandia/resources/i18n` directory. + +TODO: remove the get_relative_paths function if your stack is Python 3.12+, which added +the walk_up option in pathlib.Path.relative_to(). It would be a cleaner solution. +Listing files would become: + + python_files = [ + f.relative_to(i18n_path, walk_up=True) + for f in sorted(list(Path(src_path).rglob("*.py"))) + if not f.name.startswith("__") + ] + ui_files = [ + f.relative_to(i18n_path, walk_up=True) + for f in sorted(list(Path(src_path).rglob("*.ui"))) + ] + ts_files = [ + f.relative_to(i18n_path, walk_up=True) + for f in sorted(list(Path(src_path).rglob("*.ts"))) + ] + +So update the script accordingly by removing the function and references. +""" + +# -- Imports +from os import path +from pathlib import Path + +# -- Variables +src_path = Path("plugin_map2loop") +i18n_path = src_path.joinpath("resources/i18n") +output_file = i18n_path.joinpath("plugin_translation.pro") + +# make sure the output directory exists +i18n_path.mkdir(parents=True, exist_ok=True) + + +# -- Functions +def get_relative_paths(filepaths_list: list[Path]) -> list[str]: + """Parse a list of file paths and return a list of relative paths to the i18n + directory, escaping the backslashes to make them compatible with Qt Linguist. + + :param filepaths_list: List of file paths to parse + :type filepaths_list: list[Path] + + :return: List of relative paths to the i18n directory + :rtype: list[str] + """ + return [ + path.relpath(filepath, i18n_path).replace("\\", "/") + for filepath in filepaths_list + ] + + +# -- Run + +# Get the list of all files in directory tree at given path +python_files = [ + f for f in sorted(Path(src_path).rglob("*.py")) if not f.name.startswith("__") +] +ui_files = sorted(Path(src_path).rglob("*.ui")) +ts_files = sorted(Path(src_path).rglob("*.ts")) + + +# Generate the translation profile +forms = "FORMS =" + " \\\n".join([f"\t{f}" for f in get_relative_paths(ui_files)]) + +sources = "SOURCES =" + " \\\n".join( + [f"\t{f}" for f in get_relative_paths(python_files)] +) + +translations = "TRANSLATIONS =" + " \\\n".join( + [f"\t{f}" for f in get_relative_paths(ts_files)] +) + +# write to output file +with output_file.open("w", encoding="UTF-8") as f: + f.write( + "# This file is auto-generated by the script generate_translation_profile.py\n" + "# Do not edit this file manually\n\n" + ) + f.write(f"{forms}\n\n{sources}\n\n{translations}") diff --git a/setup.cfg b/setup.cfg index b2e1466..bb4854c 100644 --- a/setup.cfg +++ b/setup.cfg @@ -4,7 +4,7 @@ description-file = README.md [qgis-plugin-ci] plugin_path = loopstructural -project_slug = TO BE SET WITH THE GITLAB/GITHUB SLUGIFIED PROJECT NAME (IN PROJECT'S URL) +project_slug = TO BE SET WITH THE GITLAB/GITHUB SLUGIFIED PROJECT NAME (IN PROJECTS URL) github_organization_slug = TO BE SET WITH THE GITHUB USER/ORGANIZATION NAME diff --git a/tests/qgis/input/all_contacts.gpkg b/tests/qgis/input/all_contacts.gpkg new file mode 100644 index 0000000..d116bff Binary files /dev/null and b/tests/qgis/input/all_contacts.gpkg differ diff --git a/tests/qgis/input/basal_contact.gpkg b/tests/qgis/input/basal_contact.gpkg new file mode 100644 index 0000000..107a689 Binary files /dev/null and b/tests/qgis/input/basal_contact.gpkg differ diff --git a/tests/qgis/input/dtm_rp.tif b/tests/qgis/input/dtm_rp.tif new file mode 100644 index 0000000..a9ba843 Binary files /dev/null and b/tests/qgis/input/dtm_rp.tif differ diff --git a/tests/qgis/input/dtm_rp.tif.aux.xml b/tests/qgis/input/dtm_rp.tif.aux.xml new file mode 100644 index 0000000..923608b --- /dev/null +++ b/tests/qgis/input/dtm_rp.tif.aux.xml @@ -0,0 +1,11 @@ + + + + 1175 + 561.72162965205 + 297 + 107.52060579962 + 99.73 + + + diff --git a/tests/qgis/input/faults_clip.cpg b/tests/qgis/input/faults_clip.cpg new file mode 100644 index 0000000..cd89cb9 --- /dev/null +++ b/tests/qgis/input/faults_clip.cpg @@ -0,0 +1 @@ +ISO-8859-1 \ No newline at end of file diff --git a/tests/qgis/input/faults_clip.dbf b/tests/qgis/input/faults_clip.dbf new file mode 100644 index 0000000..35f7f0e Binary files /dev/null and b/tests/qgis/input/faults_clip.dbf differ diff --git a/tests/qgis/input/faults_clip.prj b/tests/qgis/input/faults_clip.prj new file mode 100644 index 0000000..51bb0e5 --- /dev/null +++ b/tests/qgis/input/faults_clip.prj @@ -0,0 +1 @@ +PROJCS["GDA94_MGA_zone_50",GEOGCS["GCS_GDA94",DATUM["D_Geocentric_Datum_of_Australia_1994",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",10000000.0],PARAMETER["Central_Meridian",117.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]] \ No newline at end of file diff --git a/tests/qgis/input/faults_clip.shp b/tests/qgis/input/faults_clip.shp new file mode 100644 index 0000000..7a0e9a2 Binary files /dev/null and b/tests/qgis/input/faults_clip.shp differ diff --git a/tests/qgis/input/faults_clip.shx b/tests/qgis/input/faults_clip.shx new file mode 100644 index 0000000..f129118 Binary files /dev/null and b/tests/qgis/input/faults_clip.shx differ diff --git a/tests/qgis/input/folds_clip.cpg b/tests/qgis/input/folds_clip.cpg new file mode 100644 index 0000000..cd89cb9 --- /dev/null +++ b/tests/qgis/input/folds_clip.cpg @@ -0,0 +1 @@ +ISO-8859-1 \ No newline at end of file diff --git a/tests/qgis/input/folds_clip.dbf b/tests/qgis/input/folds_clip.dbf new file mode 100644 index 0000000..fc04d6c Binary files /dev/null and b/tests/qgis/input/folds_clip.dbf differ diff --git a/tests/qgis/input/folds_clip.prj b/tests/qgis/input/folds_clip.prj new file mode 100644 index 0000000..51bb0e5 --- /dev/null +++ b/tests/qgis/input/folds_clip.prj @@ -0,0 +1 @@ +PROJCS["GDA94_MGA_zone_50",GEOGCS["GCS_GDA94",DATUM["D_Geocentric_Datum_of_Australia_1994",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",10000000.0],PARAMETER["Central_Meridian",117.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]] \ No newline at end of file diff --git a/tests/qgis/input/folds_clip.shp b/tests/qgis/input/folds_clip.shp new file mode 100644 index 0000000..509566c Binary files /dev/null and b/tests/qgis/input/folds_clip.shp differ diff --git a/tests/qgis/input/folds_clip.shx b/tests/qgis/input/folds_clip.shx new file mode 100644 index 0000000..981c6ce Binary files /dev/null and b/tests/qgis/input/folds_clip.shx differ diff --git a/tests/qgis/input/geol_clip_no_gaps.cpg b/tests/qgis/input/geol_clip_no_gaps.cpg new file mode 100644 index 0000000..cd89cb9 --- /dev/null +++ b/tests/qgis/input/geol_clip_no_gaps.cpg @@ -0,0 +1 @@ +ISO-8859-1 \ No newline at end of file diff --git a/tests/qgis/input/geol_clip_no_gaps.dbf b/tests/qgis/input/geol_clip_no_gaps.dbf new file mode 100644 index 0000000..832c68c Binary files /dev/null and b/tests/qgis/input/geol_clip_no_gaps.dbf differ diff --git a/tests/qgis/input/geol_clip_no_gaps.prj b/tests/qgis/input/geol_clip_no_gaps.prj new file mode 100644 index 0000000..51bb0e5 --- /dev/null +++ b/tests/qgis/input/geol_clip_no_gaps.prj @@ -0,0 +1 @@ +PROJCS["GDA94_MGA_zone_50",GEOGCS["GCS_GDA94",DATUM["D_Geocentric_Datum_of_Australia_1994",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",10000000.0],PARAMETER["Central_Meridian",117.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]] \ No newline at end of file diff --git a/tests/qgis/input/geol_clip_no_gaps.shp b/tests/qgis/input/geol_clip_no_gaps.shp new file mode 100644 index 0000000..0d7d669 Binary files /dev/null and b/tests/qgis/input/geol_clip_no_gaps.shp differ diff --git a/tests/qgis/input/geol_clip_no_gaps.shx b/tests/qgis/input/geol_clip_no_gaps.shx new file mode 100644 index 0000000..78d33f6 Binary files /dev/null and b/tests/qgis/input/geol_clip_no_gaps.shx differ diff --git a/tests/qgis/input/stratigraphic_column_testing.gpkg b/tests/qgis/input/stratigraphic_column_testing.gpkg new file mode 100644 index 0000000..bb782ef Binary files /dev/null and b/tests/qgis/input/stratigraphic_column_testing.gpkg differ diff --git a/tests/qgis/input/structure_clip.cpg b/tests/qgis/input/structure_clip.cpg new file mode 100644 index 0000000..cd89cb9 --- /dev/null +++ b/tests/qgis/input/structure_clip.cpg @@ -0,0 +1 @@ +ISO-8859-1 \ No newline at end of file diff --git a/tests/qgis/input/structure_clip.dbf b/tests/qgis/input/structure_clip.dbf new file mode 100644 index 0000000..7f15994 Binary files /dev/null and b/tests/qgis/input/structure_clip.dbf differ diff --git a/tests/qgis/input/structure_clip.prj b/tests/qgis/input/structure_clip.prj new file mode 100644 index 0000000..51bb0e5 --- /dev/null +++ b/tests/qgis/input/structure_clip.prj @@ -0,0 +1 @@ +PROJCS["GDA94_MGA_zone_50",GEOGCS["GCS_GDA94",DATUM["D_Geocentric_Datum_of_Australia_1994",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",10000000.0],PARAMETER["Central_Meridian",117.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]] \ No newline at end of file diff --git a/tests/qgis/input/structure_clip.shp b/tests/qgis/input/structure_clip.shp new file mode 100644 index 0000000..60a7773 Binary files /dev/null and b/tests/qgis/input/structure_clip.shp differ diff --git a/tests/qgis/input/structure_clip.shx b/tests/qgis/input/structure_clip.shx new file mode 100644 index 0000000..cb743f8 Binary files /dev/null and b/tests/qgis/input/structure_clip.shx differ diff --git a/tests/qgis/test_basal_contacts.py b/tests/qgis/test_basal_contacts.py new file mode 100644 index 0000000..20d0276 --- /dev/null +++ b/tests/qgis/test_basal_contacts.py @@ -0,0 +1,113 @@ +import unittest +from pathlib import Path +from qgis.core import QgsVectorLayer, QgsProcessingContext, QgsProcessingFeedback, QgsMessageLog, Qgis, QgsApplication, QgsFeature, QgsField +from qgis.PyQt.QtCore import QVariant +from qgis.testing import start_app +from loopstructural.processing.algorithms.extract_basal_contacts import BasalContactsAlgorithm +from loopstructural.processing.provider import Map2LoopProvider + +class TestBasalContacts(unittest.TestCase): + + @classmethod + def setUpClass(cls): + cls.qgs = start_app() + + cls.provider = Map2LoopProvider() + QgsApplication.processingRegistry().addProvider(cls.provider) + + def setUp(self): + self.test_dir = Path(__file__).parent + self.input_dir = self.test_dir / "input" + + self.geology_file = self.input_dir / "geol_clip_no_gaps.shp" + self.faults_file = self.input_dir / "faults_clip.shp" + self.strati_file = self.input_dir / "stratigraphic_column_testing.gpkg" + + self.assertTrue(self.geology_file.exists(), f"geology not found: {self.geology_file}") + self.assertTrue(self.strati_file.exists(), f"strati not found: {self.strati_file}") + if not self.faults_file.exists(): + QgsMessageLog.logMessage(f"faults not found: {self.faults_file}, will run test without faults", "TestBasalContacts", Qgis.Warning) + + def test_basal_contacts_extraction(self): + + geology_layer = QgsVectorLayer(str(self.geology_file), "geology", "ogr") + + self.assertTrue(geology_layer.isValid(), "geology layer should be valid") + self.assertGreater(geology_layer.featureCount(), 0, "geology layer should have features") + + faults_layer = None + if self.faults_file.exists(): + faults_layer = QgsVectorLayer(str(self.faults_file), "faults", "ogr") + self.assertTrue(faults_layer.isValid(), "faults layer should be valid") + self.assertGreater(faults_layer.featureCount(), 0, "faults layer should have features") + QgsMessageLog.logMessage(f"faults layer: {faults_layer.featureCount()} features", "TestBasalContacts", Qgis.Critical) + + QgsMessageLog.logMessage(f"geology layer: {geology_layer.featureCount()} features", "TestBasalContacts", Qgis.Critical) + + strati_table = QgsVectorLayer(str(self.strati_file), "strati", "ogr") + algorithm = BasalContactsAlgorithm() + algorithm.initAlgorithm() + + parameters = { + 'GEOLOGY': geology_layer, + 'UNIT_NAME_FIELD': 'unitname', + 'FORMATION_FIELD': 'formation', + 'FAULTS': faults_layer, + 'STRATIGRAPHIC_COLUMN': strati_table, + 'IGNORE_UNITS': [], + 'BASAL_CONTACTS': 'memory:basal_contacts', + 'ALL_CONTACTS': 'memory:all_contacts' + } + + context = QgsProcessingContext() + feedback = QgsProcessingFeedback() + + try: + QgsMessageLog.logMessage("Starting basal contacts algorithm...", "TestBasalContacts", Qgis.Critical) + + result = algorithm.processAlgorithm(parameters, context, feedback) + + QgsMessageLog.logMessage(f"Result: {result}", "TestBasalContacts", Qgis.Critical) + + self.assertIsNotNone(result, "result should not be None") + self.assertIn('BASAL_CONTACTS', result, "Result should contain BASAL_CONTACTS key") + self.assertIn('ALL_CONTACTS', result, "Result should contain ALL_CONTACTS key") + + basal_contacts_layer = context.takeResultLayer(result['BASAL_CONTACTS']) + self.assertIsNotNone(basal_contacts_layer, "basal contacts layer should not be None") + self.assertTrue(basal_contacts_layer.isValid(), "basal contacts layer should be valid") + self.assertGreater(basal_contacts_layer.featureCount(), 0, "basal contacts layer should have features") + + QgsMessageLog.logMessage(f"Generated {basal_contacts_layer.featureCount()} basal contacts", + "TestBasalContacts", Qgis.Critical) + + all_contacts_layer = context.takeResultLayer(result['ALL_CONTACTS']) + self.assertIsNotNone(all_contacts_layer, "all contacts layer should not be None") + self.assertTrue(all_contacts_layer.isValid(), "all contacts layer should be valid") + self.assertGreater(all_contacts_layer.featureCount(), 0, "all contacts layer should have features") + + QgsMessageLog.logMessage(f"Generated {all_contacts_layer.featureCount()} total contacts", + "TestBasalContacts", Qgis.Critical) + + QgsMessageLog.logMessage("Basal contacts test completed successfully!", "TestBasalContacts", Qgis.Critical) + + except Exception as e: + QgsMessageLog.logMessage(f"Basal contacts test error: {str(e)}", "TestBasalContacts", Qgis.Critical) + QgsMessageLog.logMessage(f"Error type: {type(e).__name__}", "TestBasalContacts", Qgis.Critical) + import traceback + QgsMessageLog.logMessage(f"Full traceback:\n{traceback.format_exc()}", "TestBasalContacts", Qgis.Critical) + raise + + finally: + QgsMessageLog.logMessage("=" * 50, "TestBasalContacts", Qgis.Critical) + + @classmethod + def tearDownClass(cls): + try: + registry = QgsApplication.processingRegistry() + registry.removeProvider(cls.provider) + except Exception: + pass + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/tests/qgis/test_env_var_parser.py b/tests/qgis/test_env_var_parser.py new file mode 100644 index 0000000..316b034 --- /dev/null +++ b/tests/qgis/test_env_var_parser.py @@ -0,0 +1,76 @@ +import os +import unittest + +from loopstructural.toolbelt.env_var_parser import EnvVarParser + + +class TestEnvVarParser(unittest.TestCase): + """Unit tests for EnvVarParser.""" + + def setUp(self) -> None: + """Prepare the test environment before each test.""" + self.env_backup = dict(os.environ) # Backup environment variables + + def tearDown(self) -> None: + """Restore the original environment variables after each test.""" + os.environ.clear() + os.environ.update(self.env_backup) + + def test_int_conversion(self) -> None: + """Test integer conversion from environment variable""" + os.environ["MY_INT"] = "42" + self.assertEqual(EnvVarParser.get_env_var("MY_INT", 0), 42) + + def test_float_conversion(self) -> None: + """Test float conversion from environment variable""" + os.environ["MY_FLOAT"] = "3.14" + self.assertEqual(EnvVarParser.get_env_var("MY_FLOAT", 0.0), 3.14) + + def test_bool_conversion_true(self) -> None: + """Test bool conversion from environment variable""" + for true_value in ["1", "true", "yes", "on"]: + os.environ["MY_BOOL"] = true_value + self.assertTrue(EnvVarParser.get_env_var("MY_BOOL", False)) + + def test_bool_conversion_false(self) -> None: + """Test bool conversion from environment variable""" + for false_value in ["0", "false", "no", "off"]: + os.environ["MY_BOOL"] = false_value + self.assertFalse(EnvVarParser.get_env_var("MY_BOOL", True)) + + def test_bool_invalid_defaults_to_original(self) -> None: + """Test invalid bool conversion from environment variable""" + os.environ["MY_BOOL"] = "maybe" + self.assertFalse( + EnvVarParser.get_env_var("MY_BOOL", False) + ) # Default should remain + + def test_string_conversion(self) -> None: + """Test string conversion from environment variable""" + os.environ["MY_STRING"] = "Hello, World!" + self.assertEqual( + EnvVarParser.get_env_var("MY_STRING", "default"), "Hello, World!" + ) + + def test_default_value_when_env_missing(self) -> None: + """Test default value is used when environment variable is missing""" + self.assertEqual(EnvVarParser.get_env_var("MISSING_INT", 99), 99) + + def test_invalid_int_fallback_to_default(self) -> None: + """Test default value used when the environment variable is not a valid int""" + os.environ["MY_INT"] = "not_an_int" + self.assertEqual(EnvVarParser.get_env_var("MY_INT", 10), 10) + + def test_invalid_float_fallback_to_default(self) -> None: + """Test default value used when the environment variable is not a valid float""" + os.environ["MY_FLOAT"] = "not_a_float" + self.assertEqual(EnvVarParser.get_env_var("MY_FLOAT", 1.23), 1.23) + + def test_unsupported_type(self) -> None: + """Test exception is raised when the type expected is not supported""" + os.environ["INT_LIST"] = "1,2,3,4" + self.assertRaises(TypeError, EnvVarParser.get_env_var, "INT_LIST", [1, 2, 3, 4]) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/qgis/test_processing.py b/tests/qgis/test_processing.py new file mode 100644 index 0000000..ca1a94c --- /dev/null +++ b/tests/qgis/test_processing.py @@ -0,0 +1,40 @@ +#! python3 + +""" +Usage from the repo root folder: + +.. code-block:: bash + + # for whole tests + python -m unittest tests.qgis.test_plg_processing + # for specific test + python -m unittest tests.qgis.test_plg_processing.TestPlgprocessing.test_plg_processing_structure +""" + +# PyQGIS +from qgis.core import QgsApplication +from qgis.testing import start_app, unittest + +from loopstructural.processing.provider import ( + Map2LoopProvider, +) + +provider = None + + +class TestProcessing(unittest.TestCase): + """Tests for processing algorithms.""" + + def setUp(self) -> None: + """Set up the processing tests.""" + if not QgsApplication.processingRegistry().providers(): + self.provider = Map2LoopProvider() + QgsApplication.processingRegistry().addProvider(self.provider) + self.maxDiff = None + + # Start App needed to run processing on unittest + start_app() + + def test_processing_provider(self): + """Sample test.""" + print(f"Processing provider name : {self.provider.name()}") diff --git a/tests/qgis/test_sampler_decimator.py b/tests/qgis/test_sampler_decimator.py new file mode 100644 index 0000000..c329321 --- /dev/null +++ b/tests/qgis/test_sampler_decimator.py @@ -0,0 +1,102 @@ +import unittest +from pathlib import Path +from qgis.core import QgsVectorLayer, QgsRasterLayer, QgsProcessingContext, QgsProcessingFeedback, QgsMessageLog, Qgis,QgsApplication +from qgis.testing import start_app +from loopstructural.processing.algorithms.sampler import SamplerAlgorithm +from loopstructural.processing.provider import Map2LoopProvider + +class TestSamplerDecimator(unittest.TestCase): + + @classmethod + def setUpClass(cls): + cls.qgs = start_app() + + cls.provider = Map2LoopProvider() + QgsApplication.processingRegistry().addProvider(cls.provider) + + def setUp(self): + self.test_dir = Path(__file__).parent + self.input_dir = self.test_dir / "input" + + self.geology_file = self.input_dir / "geol_clip_no_gaps.shp" + self.structure_file = self.input_dir / "structure_clip.shp" + self.dtm_file = self.input_dir / "dtm_rp.tif" + + self.assertTrue(self.geology_file.exists(), f"geology not found: {self.geology_file}") + self.assertTrue(self.structure_file.exists(), f"structure not found: {self.structure_file}") + self.assertTrue(self.dtm_file.exists(), f"dtm not found: {self.dtm_file}") + + def test_decimator_1_with_structure(self): + + geology_layer = QgsVectorLayer(str(self.geology_file), "geology", "ogr") + structure_layer = QgsVectorLayer(str(self.structure_file), "structure", "ogr") + dtm_layer = QgsRasterLayer(str(self.dtm_file), "dtm") + + self.assertTrue(geology_layer.isValid(), "geology layer should be valid") + self.assertTrue(structure_layer.isValid(), "structure layer should be valid") + self.assertTrue(dtm_layer.isValid(), "dtm layer should be valid") + self.assertGreater(geology_layer.featureCount(), 0, "geology layer should have features") + self.assertGreater(structure_layer.featureCount(), 0, "structure layer should have features") + + QgsMessageLog.logMessage(f"geology layer valid: {geology_layer.isValid()}", "TestDecimator", Qgis.Critical) + QgsMessageLog.logMessage(f"structure layer valid: {structure_layer.isValid()}", "TestDecimator", Qgis.Critical) + QgsMessageLog.logMessage(f"dtm layer valid: {dtm_layer.isValid()}", "TestDecimator", Qgis.Critical) + QgsMessageLog.logMessage(f"dtm source: {dtm_layer.source()}", "TestDecimator", Qgis.Critical) + + QgsMessageLog.logMessage(f"geology layer: {geology_layer.featureCount()} features", "TestDecimator", Qgis.Critical) + QgsMessageLog.logMessage(f"structure layer: {structure_layer.featureCount()} features", "TestDecimator", Qgis.Critical) + QgsMessageLog.logMessage(f"spatial data- structure layer", "TestDecimator", Qgis.Critical) + QgsMessageLog.logMessage(f"sampler type: Decimator", "TestDecimator", Qgis.Critical) + QgsMessageLog.logMessage(f"decimation: 1", "TestDecimator", Qgis.Critical) + QgsMessageLog.logMessage(f"dtm: {self.dtm_file.name}", "TestDecimator", Qgis.Critical) + + algorithm = SamplerAlgorithm() + algorithm.initAlgorithm() + + parameters = { + 'DTM': dtm_layer, + 'GEOLOGY': geology_layer, + 'SPATIAL_DATA': structure_layer, + 'SAMPLER_TYPE': 0, + 'DECIMATION': 1, + 'SPACING': 200.0, + 'SAMPLED_CONTACTS': 'memory:decimated_points' + } + + context = QgsProcessingContext() + feedback = QgsProcessingFeedback() + + + try: + QgsMessageLog.logMessage("Starting decimator sampler algorithm...", "TestDecimator", Qgis.Critical) + + result = algorithm.processAlgorithm(parameters, context, feedback) + + QgsMessageLog.logMessage(f"Result: {result}", "TestDecimator", Qgis.Critical) + + self.assertIsNotNone(result, "result should not be None") + self.assertIn('SAMPLED_CONTACTS', result, "Result should contain SAMPLED_CONTACTS key") + + QgsMessageLog.logMessage("Decimator sampler test completed successfully!", "TestDecimator", Qgis.Critical) + + except Exception as e: + QgsMessageLog.logMessage(f"Decimator sampler test error: {str(e)}", "TestDecimator", Qgis.Critical) + QgsMessageLog.logMessage(f"Error type: {type(e).__name__}", "TestDecimator", Qgis.Critical) + + import traceback + QgsMessageLog.logMessage(f"Full traceback:\n{traceback.format_exc()}", "TestDecimator", Qgis.Critical) + raise + + finally: + QgsMessageLog.logMessage("=" * 50, "TestDecimator", Qgis.Critical) + + @classmethod + def tearDownClass(cls): + try: + registry = QgsApplication.processingRegistry() + registry.removeProvider(cls.provider) + except Exception: + pass + +if __name__ == '__main__': + unittest.main() diff --git a/tests/qgis/test_sampler_spacing.py b/tests/qgis/test_sampler_spacing.py new file mode 100644 index 0000000..e2f285d --- /dev/null +++ b/tests/qgis/test_sampler_spacing.py @@ -0,0 +1,84 @@ +import unittest +from pathlib import Path +from qgis.core import QgsVectorLayer, QgsProcessingContext, QgsProcessingFeedback, QgsMessageLog, Qgis, QgsApplication +from qgis.testing import start_app +from loopstructural.processing.algorithms.sampler import SamplerAlgorithm +from loopstructural.processing.provider import Map2LoopProvider + +class TestSamplerSpacing(unittest.TestCase): + + @classmethod + def setUpClass(cls): + cls.qgs = start_app() + + cls.provider = Map2LoopProvider() + QgsApplication.processingRegistry().addProvider(cls.provider) + + def setUp(self): + self.test_dir = Path(__file__).parent + self.input_dir = self.test_dir / "input" + + self.geology_file = self.input_dir / "geol_clip_no_gaps.shp" + + self.assertTrue(self.geology_file.exists(), f"geology not found: {self.geology_file}") + + def test_spacing_50_with_geology(self): + + geology_layer = QgsVectorLayer(str(self.geology_file), "geology", "ogr") + + self.assertTrue(geology_layer.isValid(), "geology layer should be valid") + self.assertGreater(geology_layer.featureCount(), 0, "geology layer should have features") + + QgsMessageLog.logMessage(f"geology layer: {geology_layer.featureCount()} features", "TestSampler", Qgis.Critical) + QgsMessageLog.logMessage(f"spatial data- geology layer", "TestSampler", Qgis.Critical) + QgsMessageLog.logMessage(f"sampler type: Spacing", "TestSampler", Qgis.Critical) + QgsMessageLog.logMessage(f"spacing: 50", "TestSampler", Qgis.Critical) + + algorithm = SamplerAlgorithm() + algorithm.initAlgorithm() + + parameters = { + 'DTM': None, + 'GEOLOGY': None, + 'SPATIAL_DATA': geology_layer, + 'SAMPLER_TYPE': 1, + 'DECIMATION': 1, + 'SPACING': 50.0, + 'SAMPLED_CONTACTS': 'memory:sampled_points' + } + + context = QgsProcessingContext() + feedback = QgsProcessingFeedback() + + try: + QgsMessageLog.logMessage("Starting spacing sampler algorithm...", "TestSampler", Qgis.Critical) + + result = algorithm.processAlgorithm(parameters, context, feedback) + + QgsMessageLog.logMessage(f"Result: {result}", "TestSampler", Qgis.Critical) + + self.assertIsNotNone(result, "result should not be None") + self.assertIn('SAMPLED_CONTACTS', result, "Result should contain SAMPLED_CONTACTS key") + + QgsMessageLog.logMessage("Spacing sampler test completed successfully!", "TestSampler", Qgis.Critical) + + except Exception as e: + QgsMessageLog.logMessage(f"Spacing sampler test error: {str(e)}", "TestSampler", Qgis.Critical) + QgsMessageLog.logMessage(f"Error type: {type(e).__name__}", "TestSampler", Qgis.Critical) + import traceback + QgsMessageLog.logMessage(f"Full traceback:\n{traceback.format_exc()}", "TestSampler", Qgis.Critical) + raise + + finally: + QgsMessageLog.logMessage("=" * 50, "TestSampler", Qgis.Critical) + + @classmethod + def tearDownClass(cls): + try: + registry = QgsApplication.processingRegistry() + registry.removeProvider(cls.provider) + except Exception: + pass + +if __name__ == '__main__': + unittest.main()