diff --git a/.github/workflows/linter.yml b/.github/workflows/linter.yml index 57c2289..a725d82 100644 --- a/.github/workflows/linter.yml +++ b/.github/workflows/linter.yml @@ -7,13 +7,13 @@ on: pull_request: branches: - - master + - main paths: - '**.py' workflow_dispatch: env: - PROJECT_FOLDER: "map2loop" + PROJECT_FOLDER: "m2l" PYTHON_VERSION: 3.9 permissions: contents: write @@ -55,5 +55,5 @@ jobs: token: ${{ secrets.GITHUB_TOKEN }} title: "style: auto format fixes" body: "This PR applies style fixes by black and ruff." - base: master + base: main branch: lint/style-fixes-${{ github.run_id }} diff --git a/.github/workflows/tester.yml b/.github/workflows/tester.yml index 4543564..cf13a1b 100644 --- a/.github/workflows/tester.yml +++ b/.github/workflows/tester.yml @@ -2,16 +2,12 @@ name: "🎳 Tester" on: push: - branches: - - main paths: - '**.py' - .github/workflows/tester.yml - requirements/testing.txt pull_request: - branches: - - main paths: - '**.py' - .github/workflows/tester.yml @@ -45,55 +41,62 @@ jobs: - name: Run Unit tests run: pytest -p no:qgis tests/unit/ - # test-qgis: - # runs-on: ubuntu-latest - - # container: - # image: qgis/qgis:3.4 - # env: - # CI: true - # DISPLAY: ":1" - # MUTE_LOGS: true - # NO_MODALS: 1 - # PYTHONPATH: "/usr/share/qgis/python/plugins:/usr/share/qgis/python:." - # QT_QPA_PLATFORM: "offscreen" - # WITH_PYTHON_PEP: false - # # be careful, things have changed since QGIS 3.40. So if you are using this setup - # # with a QGIS version older than 3.40, you may need to change the way you set up the container - # volumes: - # # Mount the X11 socket to allow GUI applications to run - # - /tmp/.X11-unix:/tmp/.X11-unix - # # Mount the workspace directory to the container - # - ${{ github.workspace }}:/home/root/ - - # steps: - # - name: Get source code - # uses: actions/checkout@v4 - - # - name: Print QGIS version - # run: qgis --version - - # # Uncomment if you need to run a script to set up the plugin in QGIS docker image < 3.40 - # # - name: Setup plugin - # # run: qgis_setup.sh ${{ env.PROJECT_FOLDER }} - - # - name: Install Python requirements - # run: | - # apt update && apt install -y python3-pip python3-venv pipx - # # Create a virtual environment - # cd /home/root/ - # pipx run qgis-venv-creator --venv-name ".venv" - # # Activate the virtual environment - # . .venv/bin/activate - # # Install the requirements - # python3 -m pip install -U -r requirements/testing.txt - - # - name: Run Unit tests - # run: | - # cd /home/root/ - # # Activate the virtual environment - # . .venv/bin/activate - # # Run the tests - # # xvfb-run is used to run the tests in a virtual framebuffer - # # This is necessary because QGIS requires a display to run - # xvfb-run python3 -m pytest tests/qgis --junitxml=junit/test-results-qgis.xml --cov-report=xml:coverage-reports/coverage-qgis.xml + test-qgis: + runs-on: ubuntu-latest + + container: + image: qgis/qgis:latest + env: + CI: true + DISPLAY: ":1" + MUTE_LOGS: true + NO_MODALS: 1 + PYTHONPATH: "/usr/share/qgis/python/plugins:/usr/share/qgis/python:." + QT_QPA_PLATFORM: "offscreen" + WITH_PYTHON_PEP: false + # be careful, things have changed since QGIS 3.40. So if you are using this setup + # with a QGIS version older than 3.40, you may need to change the way you set up the container + volumes: + # Mount the X11 socket to allow GUI applications to run + - /tmp/.X11-unix:/tmp/.X11-unix + # Mount the workspace directory to the container + - ${{ github.workspace }}:/home/root/ + + steps: + - name: Get source code + uses: actions/checkout@v4 + + - name: Print QGIS version + run: qgis --version + + # Uncomment if you need to run a script to set up the plugin in QGIS docker image < 3.40 + # - name: Setup plugin + # run: qgis_setup.sh ${{ env.PROJECT_FOLDER }} + + - name: Install Python requirements + run: | + apt update && apt install -y python3-pip python3-venv pipx + # Create a virtual environment + cd /home/root/ + pipx run qgis-venv-creator --venv-name ".venv" + # Activate the virtual environment + . .venv/bin/activate + # Install the requirements + python3 -m pip install -U -r requirements/testing.txt + python3 -m pip install git+https://github.com/Loop3D/map2loop.git@noelle/contact_extractor + + - name: verify input data + run: | + cd /home/root/ + . .venv/bin/activate + ls -la tests/qgis/input/ || echo "Input directory not found" + + - name: Run Unit tests + run: | + cd /home/root/ + # Activate the virtual environment + . .venv/bin/activate + # Run the tests + # xvfb-run is used to run the tests in a virtual framebuffer + # This is necessary because QGIS requires a display to run + xvfb-run python3 -m pytest tests/qgis --junitxml=junit/test-results-qgis.xml --cov-report=xml:coverage-reports/coverage-qgis.xml diff --git a/docs/conf.py b/docs/conf.py index 6d7892e..d7cc5df 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -15,7 +15,7 @@ import sphinx_rtd_theme # noqa: F401 theme of Read the Docs # Package -from map2loop import __about__ +from m2l import __about__ # -- Build environment ----------------------------------------------------- on_rtd = environ.get("READTHEDOCS", None) == "True" diff --git a/map2loop/__about__.py b/m2l/__about__.py similarity index 100% rename from map2loop/__about__.py rename to m2l/__about__.py diff --git a/map2loop/__init__.py b/m2l/__init__.py similarity index 100% rename from map2loop/__init__.py rename to m2l/__init__.py diff --git a/map2loop/gui/__init__.py b/m2l/gui/__init__.py similarity index 100% rename from map2loop/gui/__init__.py rename to m2l/gui/__init__.py diff --git a/map2loop/gui/dlg_settings.py b/m2l/gui/dlg_settings.py similarity index 96% rename from map2loop/gui/dlg_settings.py rename to m2l/gui/dlg_settings.py index bf6c8b1..351a008 100644 --- a/map2loop/gui/dlg_settings.py +++ b/m2l/gui/dlg_settings.py @@ -18,15 +18,15 @@ from qgis.PyQt.QtGui import QDesktopServices, QIcon # project -from map2loop.__about__ import ( +from m2l.__about__ import ( __icon_path__, __title__, __uri_homepage__, __uri_tracker__, __version__, ) -from map2loop.toolbelt import PlgLogger, PlgOptionsManager -from map2loop.toolbelt.preferences import PlgSettingsStructure +from m2l.toolbelt import PlgLogger, PlgOptionsManager +from m2l.toolbelt.preferences import PlgSettingsStructure # ############################################################################ # ########## Globals ############### diff --git a/map2loop/gui/dlg_settings.ui b/m2l/gui/dlg_settings.ui similarity index 100% rename from map2loop/gui/dlg_settings.ui rename to m2l/gui/dlg_settings.ui diff --git a/m2l/main/vectorLayerWrapper.py b/m2l/main/vectorLayerWrapper.py new file mode 100644 index 0000000..4476e6a --- /dev/null +++ b/m2l/main/vectorLayerWrapper.py @@ -0,0 +1,727 @@ +# PyQGIS / PyQt imports +from osgeo import gdal +from qgis.core import ( + QgsRaster, + QgsFields, + QgsField, + QgsFeature, + QgsGeometry, + QgsWkbTypes, + QgsCoordinateReferenceSystem, + QgsFeatureSink, + QgsProcessingException, + QgsPoint, + QgsPointXY, + QgsProject, + QgsCoordinateTransform, + QgsRasterLayer + ) + +from qgis.PyQt.QtCore import QVariant, QDateTime +from qgis import processing +from shapely.geometry import Point, MultiPoint, LineString, MultiLineString, Polygon, MultiPolygon +from shapely.wkb import loads as wkb_loads +import pandas as pd +import geopandas as gpd +import numpy as np +import tempfile +import os + +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: + if layer is None: + return None + features = layer.getFeatures() + fields = layer.fields() + data = {'geometry': []} + for f in fields: + data[f.name()] = [] + for feature in features: + geom = feature.geometry() + if geom.isEmpty(): + continue + data['geometry'].append(geom) + for f in fields: + 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] + + for f in feats: + geom = f.geometry() + pts = vertices_from_geometry(geom) + + 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 + ---------- + alg : QgsProcessingAlgorithm (self) + gdf : geopandas.GeoDataFrame + parameters : dict (from processAlgorithm) + context : QgsProcessingContext + output_key : str (e.g. self.OUTPUT) + feedback : QgsProcessingFeedback | None + + Returns + ------- + str : dest_id to return from processAlgorithm, e.g. { output_key: dest_id } + """ + + if feedback is None: + class _Dummy: + def pushInfo(self, *a, **k): pass + def reportError(self, *a, **k): pass + def setProgress(self, *a, **k): pass + def isCanceled(self): return False + feedback = _Dummy() + + if geodataframe is None: + raise ValueError("GeoDataFrame is None") + if geodataframe.empty: + feedback.pushInfo("Input GeoDataFrame is empty; creating empty output layer.") + + # --- infer WKB type (family, Multi, Z) + def _infer_wkb(series): + base = None + any_multi = False + has_z = False + for geom in series: + if geom is None: continue + if getattr(geom, "is_empty", False): continue + # multi? + if isinstance(geom, (MultiPoint, MultiLineString, MultiPolygon)): + any_multi = True + g0 = next(iter(getattr(geom, "geoms", [])), None) + gt = getattr(g0, "geom_type", None) or None + else: + gt = getattr(geom, "geom_type", None) + + # base family + if gt in ("Point", "LineString", "Polygon"): + base = gt + # z? + try: + has_z = has_z or bool(getattr(geom, "has_z", False)) + except Exception: + pass + if base: + break + + if base is None: + # default safely to LineString if everything is empty; adjust if you prefer Point/Polygon + base = "LineString" + + fam = { + "Point": QgsWkbTypes.Point, + "LineString": QgsWkbTypes.LineString, + "Polygon": QgsWkbTypes.Polygon, + }[base] + + if any_multi: + fam = QgsWkbTypes.multiType(fam) + if has_z: + fam = QgsWkbTypes.addZ(fam) + return fam + + wkb_type = _infer_wkb(geodataframe.geometry) + + # --- build CRS from gdf.crs + crs = QgsCoordinateReferenceSystem() + if geodataframe.crs is not None: + try: + crs = QgsCoordinateReferenceSystem.fromWkt(geodataframe.crs.to_wkt()) + except Exception: + try: + epsg = geodataframe.crs.to_epsg() + if epsg: + crs = QgsCoordinateReferenceSystem.fromEpsgId(int(epsg)) + except Exception: + pass + + # --- build QGIS fields from pandas dtypes + fields = QgsFields() + non_geom_cols = [c for c in geodataframe.columns if c != geodataframe.geometry.name] + + def _qvariant_type(dtype) -> QVariant.Type: + if pd.api.types.is_integer_dtype(dtype): + return QVariant.Int + if pd.api.types.is_float_dtype(dtype): + return QVariant.Double + if pd.api.types.is_bool_dtype(dtype): + return QVariant.Bool + if pd.api.types.is_datetime64_any_dtype(dtype): + return QVariant.DateTime + return QVariant.String + + for col in non_geom_cols: + fields.append(QgsField(str(col), _qvariant_type(geodataframe[col].dtype))) + + # --- create sink + sink, dest_id = qgs_algorithm.parameterAsSink( + parameters, + output_key, + context, + fields, + wkb_type, + crs, + ) + if sink is None: + raise QgsProcessingException("Could not create output sink") + + # --- write features + total = len(geodataframe.index) + is_multi_sink = QgsWkbTypes.isMultiType(wkb_type) + + for i, (_, row) in enumerate(geodataframe.iterrows()): + if feedback.isCanceled(): + break + + geom = row[geodataframe.geometry.name] + if geom is None or getattr(geom, "is_empty", False): + continue + + # promote single β†’ multi if needed + if is_multi_sink: + if isinstance(geom, Point): + geom = MultiPoint([geom]) + elif isinstance(geom, LineString): + geom = MultiLineString([geom]) + elif isinstance(geom, Polygon): + geom = MultiPolygon([geom]) + + f = QgsFeature(fields) + + # attributes in declared order + attrs = [] + for col in non_geom_cols: + val = row[col] + if isinstance(val, np.generic): + try: + val = val.item() + except Exception: + pass + if pd.api.types.is_datetime64_any_dtype(geodataframe[col].dtype): + if pd.isna(val): + val = None + else: + val = QDateTime(val.to_pydatetime()) + attrs.append(val) + f.setAttributes(attrs) + + # geometry (shapely β†’ QGIS) + try: + f.setGeometry(QgsGeometry.fromWkb(geom.wkb)) + except Exception: + f.setGeometry(QgsGeometry.fromWkt(geom.wkt)) + + sink.addFeature(f, QgsFeatureSink.FastInsert) + + if total: + feedback.setProgress(int(100.0 * (i + 1) / total)) + + return dest_id + + +# ---------- helpers ---------- + +def _qvariant_type_from_dtype(dtype) -> QVariant.Type: + """Map a pandas dtype to a QVariant type.""" + import numpy as np + if np.issubdtype(dtype, np.integer): + # prefer 64-bit when detected + try: + return QVariant.LongLong + except AttributeError: + return QVariant.Int + if np.issubdtype(dtype, np.floating): + return QVariant.Double + if np.issubdtype(dtype, np.bool_): + return QVariant.Bool + # datetimes + try: + import pandas as pd + if pd.api.types.is_datetime64_any_dtype(dtype): + return QVariant.DateTime + if pd.api.types.is_datetime64_ns_dtype(dtype): + return QVariant.DateTime + if pd.api.types.is_datetime64_dtype(dtype): + return QVariant.DateTime + if pd.api.types.is_timedelta64_dtype(dtype): + # store as string "HH:MM:SS" fallback + return QVariant.String + except Exception: + pass + # default to string + return QVariant.String + + +def _fields_from_dataframe(df, drop_cols=None) -> QgsFields: + """Build QgsFields from DataFrame dtypes.""" + drop_cols = set(drop_cols or []) + fields = QgsFields() + for name, dtype in df.dtypes.items(): + if name in drop_cols: + continue + vtype = _qvariant_type_from_dtype(dtype) + fields.append(QgsField(name, vtype)) + return fields + + +# ---------- main function you'll call inside processAlgorithm ---------- + +def dataframeToQgsLayer( + df, + x_col: str, + y_col: str, + *, + crs: QgsCoordinateReferenceSystem, + algorithm, # `self` inside a QgsProcessingAlgorithm + parameters: dict, + context, + feedback, + sink_param_name: str = "OUTPUT", + z_col: str = None, + m_col: str = None, + include_coords_in_attrs: bool = False, +): + """ + Write a pandas DataFrame to a point feature sink (QgsProcessingParameterFeatureSink). + + Params + ------ + df : pandas.DataFrame Data with coordinate columns. + x_col, y_col : str Column names for X/Easting/Longitude and Y/Northing/Latitude. + crs : QgsCoordinateReferenceSystem CRS of the coordinates (e.g., QgsCoordinateReferenceSystem('EPSG:4326')). + algorithm : QgsProcessingAlgorithm Use `self` from inside processAlgorithm. + parameters, context, feedback Standard Processing plumbing. + sink_param_name : str Name of your sink output parameter (default "OUTPUT"). + z_col, m_col : str | None Optional Z and M columns for 3D/M points. + include_coords_in_attrs : bool If False, x/y/z/m are not written as attributes. + + Returns + ------- + (sink, sink_id) The created sink and its ID. Also returns feature count via feedback. + """ + import pandas as pd + if not isinstance(df, pd.DataFrame): + raise TypeError("df must be a pandas.DataFrame") + + # Make a working copy; optionally drop coordinate columns from attributes + attr_df = df.copy() + drop_cols = [] + for col in [x_col, y_col, z_col, m_col]: + if col and not include_coords_in_attrs: + drop_cols.append(col) + + fields = _fields_from_dataframe(attr_df, drop_cols=drop_cols) + + # Geometry type (2D/3D/M) + has_z = z_col is not None and z_col in df.columns + has_m = m_col is not None and m_col in df.columns + if has_z and has_m: + wkb = QgsWkbTypes.PointZM + elif has_z: + wkb = QgsWkbTypes.PointZ + elif has_m: + wkb = QgsWkbTypes.PointM + else: + wkb = QgsWkbTypes.Point + + # Create the sink + sink, sink_id = algorithm.parameterAsSink( + parameters, + sink_param_name, + context, + fields, + wkb, + crs + ) + if sink is None: + raise QgsProcessingException("Could not create feature sink. Check output parameter and inputs.") + + total = len(df) + feedback.pushInfo(f"Writing {total} features…") + + # Precompute attribute column order + attr_columns = [f.name() for f in fields] + + # Iterate rows and write features + for i, (_idx, row) in enumerate(df.iterrows(), start=1): + if feedback.isCanceled(): + break + + # Build point geometry + x = row[x_col] + y = row[y_col] + + # skip rows with missing coords + if pd.isna(x) or pd.isna(y): + continue + + if has_z and not pd.isna(row[z_col]) and has_m and not pd.isna(row[m_col]): + pt = QgsPoint(float(x), float(y), float(row[z_col]), float(row[m_col])) + elif has_z and not pd.isna(row[z_col]): + pt = QgsPoint(float(x), float(y), float(row[z_col])) + elif has_m and not pd.isna(row[m_col]): + # PointM constructor: setZValue not needed; M is the 4th ordinate + pt = QgsPoint(float(x), float(y)) + pt.setM(float(row[m_col])) + else: + pt = QgsPointXY(float(x), float(y)) + + feat = QgsFeature(fields) + feat.setGeometry(QgsGeometry.fromPoint(pt) if isinstance(pt, QgsPoint) else QgsGeometry.fromPointXY(pt)) + + # Attributes in the same order as fields + attrs = [] + for col in attr_columns: + val = row[col] if col in row else None + # Pandas NaN -> None + if pd.isna(val): + val = None + # Convert numpy types to Python scalars to avoid QVariant issues + try: + import numpy as np + if isinstance(val, (np.generic,)): + val = val.item() + except Exception: + pass + # Convert pandas Timestamp to Python datetime + if hasattr(val, "to_pydatetime"): + try: + val = val.to_pydatetime() + except Exception: + val = str(val) + attrs.append(val) + feat.setAttributes(attrs) + + sink.addFeature(feat, QgsFeature.FastInsert) + + if i % 1000 == 0: + feedback.setProgress(int(100.0 * i / max(total, 1))) + + feedback.pushInfo("Done.") + feedback.setProgress(100) + return sink, sink_id + +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/map2loop/metadata.txt b/m2l/metadata.txt similarity index 100% rename from map2loop/metadata.txt rename to m2l/metadata.txt diff --git a/map2loop/plugin_main.py b/m2l/plugin_main.py similarity index 96% rename from map2loop/plugin_main.py rename to m2l/plugin_main.py index f803c9a..a286f1e 100644 --- a/map2loop/plugin_main.py +++ b/m2l/plugin_main.py @@ -15,17 +15,17 @@ from qgis.PyQt.QtWidgets import QAction # project -from map2loop.__about__ import ( +from m2l.__about__ import ( DIR_PLUGIN_ROOT, __icon_path__, __title__, __uri_homepage__, ) -from map2loop.gui.dlg_settings import PlgOptionsFactory -from map2loop.processing import ( +from m2l.gui.dlg_settings import PlgOptionsFactory +from m2l.processing import ( Map2LoopProvider, ) -from map2loop.toolbelt import PlgLogger +from m2l.toolbelt import PlgLogger # ############################################################################ # ########## Classes ############### diff --git a/map2loop/processing/__init__.py b/m2l/processing/__init__.py similarity index 100% rename from map2loop/processing/__init__.py rename to m2l/processing/__init__.py diff --git a/m2l/processing/algorithms/__init__.py b/m2l/processing/algorithms/__init__.py new file mode 100644 index 0000000..08d76a0 --- /dev/null +++ b/m2l/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/m2l/processing/algorithms/extract_basal_contacts.py b/m2l/processing/algorithms/extract_basal_contacts.py new file mode 100644 index 0000000..4d3ba5f --- /dev/null +++ b/m2l/processing/algorithms/extract_basal_contacts.py @@ -0,0 +1,205 @@ +""" +*************************************************************************** +* * +* 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( + QgsProcessingParameterField( + 'FORMATION_FIELD', + 'Formation Field', + parentLayerParameterName=self.INPUT_GEOLOGY, + type=QgsProcessingParameterField.String, + defaultValue='formation', + optional=True + ) + ) + + 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/m2l/processing/algorithms/sampler.py b/m2l/processing/algorithms/sampler.py new file mode 100644 index 0000000..fdc9856 --- /dev/null +++ b/m2l/processing/algorithms/sampler.py @@ -0,0 +1,240 @@ +""" +*************************************************************************** +* * +* 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 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", "Spacing"], + 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", + QgsProcessingParameterNumber.Integer, + defaultValue=1, + optional=True, + ) + ) + + self.addParameter( + QgsProcessingParameterNumber( + self.INPUT_SPACING, + "SPACING", + 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() + if samples is not None and not samples.empty: + for column_name in samples.columns: + dtype = samples[column_name].dtype + dtype_str = str(dtype) + + if dtype_str in ['float16', 'float32', 'float64']: + field_type = QVariant.Double + elif dtype_str in ['int8', 'int16', 'int32', 'int64']: + field_type = QVariant.Int + else: + field_type = QVariant.String + + fields.append(QgsField(column_name, field_type)) + + 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']))) + + attributes = [] + for column_name in samples.columns: + value = row.get(column_name) + dtype = samples[column_name].dtype + + if pd.isna(value): + attributes.append(None) + elif dtype in ['float16', 'float32', 'float64']: + attributes.append(float(value)) + elif dtype in ['int8', 'int16', 'int32', 'int64']: + attributes.append(int(value)) + else: + attributes.append(str(value)) + + feature.setAttributes(attributes) + + 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/m2l/processing/algorithms/sorter.py b/m2l/processing/algorithms/sorter.py new file mode 100644 index 0000000..ec7d086 --- /dev/null +++ b/m2l/processing/algorithms/sorter.py @@ -0,0 +1,427 @@ +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=False + ) + ) + + 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=True, + ) + ) + + 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] + sorter_name = list(SORTER_LIST.keys())[algo_index] + requires_contacts = sorter_cls in [SorterAlpha, SorterMaximiseContacts, SorterObservationProjections] + 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) + + if requires_contacts and not contacts_layer or not contacts_layer.isValid(): + raise QgsProcessingException(f"{sorter_name} requires a contacts layer") + + 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") + + unique_units = {} + units_records = [] + for f in layer.getFeatures(): + unit_name = f[unit_name_field] + if unit_name not in unique_units: + unique_units[unit_name] = len(unique_units) + units_records.append( + dict( + layerId=len(units_records), + 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/m2l/processing/algorithms/thickness_calculator.py b/m2l/processing/algorithms/thickness_calculator.py new file mode 100644 index 0000000..8d67dc5 --- /dev/null +++ b/m2l/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/m2l/processing/algorithms/user_defined_sorter.py b/m2l/processing/algorithms/user_defined_sorter.py new file mode 100644 index 0000000..23857a3 --- /dev/null +++ b/m2l/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/map2loop/processing/provider.py b/m2l/processing/provider.py similarity index 84% rename from map2loop/processing/provider.py rename to m2l/processing/provider.py index 6e0add8..9616d9b 100644 --- a/map2loop/processing/provider.py +++ b/m2l/processing/provider.py @@ -10,12 +10,20 @@ from qgis.PyQt.QtGui import QIcon # project -from map2loop.__about__ import ( +from m2l.__about__ import ( __icon_path__, __title__, __version__, ) +from .algorithms import ( + BasalContactsAlgorithm, + StratigraphySorterAlgorithm, + UserDefinedStratigraphyAlgorithm, + ThicknessCalculatorAlgorithm, + SamplerAlgorithm +) + # ############################################################################ # ########## Classes ############### # ################################## @@ -26,7 +34,11 @@ class Map2LoopProvider(QgsProcessingProvider): def loadAlgorithms(self): """Loads all algorithms belonging to this provider.""" - pass + 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, \ diff --git a/map2loop/resources/i18n/plugin_map2loop_en.ts b/m2l/resources/i18n/plugin_map2loop_en.ts similarity index 100% rename from map2loop/resources/i18n/plugin_map2loop_en.ts rename to m2l/resources/i18n/plugin_map2loop_en.ts diff --git a/map2loop/resources/i18n/plugin_translation.pro b/m2l/resources/i18n/plugin_translation.pro similarity index 100% rename from map2loop/resources/i18n/plugin_translation.pro rename to m2l/resources/i18n/plugin_translation.pro diff --git a/map2loop/resources/images/default_icon.png b/m2l/resources/images/default_icon.png similarity index 100% rename from map2loop/resources/images/default_icon.png rename to m2l/resources/images/default_icon.png diff --git a/map2loop/toolbelt/__init__.py b/m2l/toolbelt/__init__.py similarity index 100% rename from map2loop/toolbelt/__init__.py rename to m2l/toolbelt/__init__.py diff --git a/map2loop/toolbelt/env_var_parser.py b/m2l/toolbelt/env_var_parser.py similarity index 100% rename from map2loop/toolbelt/env_var_parser.py rename to m2l/toolbelt/env_var_parser.py diff --git a/map2loop/toolbelt/log_handler.py b/m2l/toolbelt/log_handler.py similarity index 98% rename from map2loop/toolbelt/log_handler.py rename to m2l/toolbelt/log_handler.py index 222fa3f..284365e 100644 --- a/map2loop/toolbelt/log_handler.py +++ b/m2l/toolbelt/log_handler.py @@ -12,8 +12,8 @@ from qgis.utils import iface # project package -import map2loop.toolbelt.preferences as plg_prefs_hdlr -from map2loop.__about__ import __title__ +import m2l.toolbelt.preferences as plg_prefs_hdlr +from m2l.__about__ import __title__ # ############################################################################ # ########## Classes ############### diff --git a/map2loop/toolbelt/preferences.py b/m2l/toolbelt/preferences.py similarity index 97% rename from map2loop/toolbelt/preferences.py rename to m2l/toolbelt/preferences.py index 85986b1..8742313 100644 --- a/map2loop/toolbelt/preferences.py +++ b/m2l/toolbelt/preferences.py @@ -11,9 +11,9 @@ from qgis.core import QgsSettings # package -import map2loop.toolbelt.log_handler as log_hdlr -from map2loop.__about__ import __title__, __version__ -from map2loop.toolbelt.env_var_parser import EnvVarParser +import m2l.toolbelt.log_handler as log_hdlr +from m2l.__about__ import __title__, __version__ +from m2l.toolbelt.env_var_parser import EnvVarParser # ############################################################################ # ########## Classes ############### diff --git a/requirements/testing.txt b/requirements/testing.txt index 1940035..d238575 100644 --- a/requirements/testing.txt +++ b/requirements/testing.txt @@ -3,3 +3,5 @@ pytest-cov>=4 packaging>=23 +shapely +geopandas \ No newline at end of file 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..bb90de8 --- /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 m2l.processing.algorithms.extract_basal_contacts import BasalContactsAlgorithm +from m2l.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 index ce4b73d..b968008 100644 --- a/tests/qgis/test_env_var_parser.py +++ b/tests/qgis/test_env_var_parser.py @@ -1,7 +1,7 @@ import os import unittest -from map2loop.toolbelt.env_var_parser import EnvVarParser +from m2l.toolbelt.env_var_parser import EnvVarParser class TestEnvVarParser(unittest.TestCase): diff --git a/tests/qgis/test_plg_preferences.py b/tests/qgis/test_plg_preferences.py index 23daa19..8ac07a3 100644 --- a/tests/qgis/test_plg_preferences.py +++ b/tests/qgis/test_plg_preferences.py @@ -18,8 +18,8 @@ from qgis.testing import unittest # project -from map2loop.__about__ import __version__ -from map2loop.toolbelt.preferences import ( +from m2l.__about__ import __version__ +from m2l.toolbelt.preferences import ( PREFIX_ENV_VARIABLE, PlgOptionsManager, PlgSettingsStructure, diff --git a/tests/qgis/test_processing.py b/tests/qgis/test_processing.py index b52b3de..6f5c719 100644 --- a/tests/qgis/test_processing.py +++ b/tests/qgis/test_processing.py @@ -15,7 +15,7 @@ from qgis.core import QgsApplication from qgis.testing import start_app, unittest -from map2loop.processing.provider import ( +from m2l.processing.provider import ( Map2LoopProvider, ) diff --git a/tests/qgis/test_sampler_decimator.py b/tests/qgis/test_sampler_decimator.py new file mode 100644 index 0000000..bb1864c --- /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 m2l.processing.algorithms.sampler import SamplerAlgorithm +from m2l.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..a49042f --- /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 m2l.processing.algorithms.sampler import SamplerAlgorithm +from m2l.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() diff --git a/tests/unit/test_plg_metadata.py b/tests/unit/test_plg_metadata.py index 2baf6ed..e6373a6 100644 --- a/tests/unit/test_plg_metadata.py +++ b/tests/unit/test_plg_metadata.py @@ -18,7 +18,7 @@ from packaging.version import parse # project -from map2loop import __about__ +from m2l import __about__ # ############################################################################ # ########## Classes #############