diff --git a/loopstructural/gui/modelling/geological_model_tab/geological_model_tab.py b/loopstructural/gui/modelling/geological_model_tab/geological_model_tab.py index 377a8d7..7cc00a5 100644 --- a/loopstructural/gui/modelling/geological_model_tab/geological_model_tab.py +++ b/loopstructural/gui/modelling/geological_model_tab/geological_model_tab.py @@ -150,6 +150,26 @@ def initialize_model(self): "Please set the bounding box before initializing the model.", ) return + + # Validate model CRS + if not self.data_manager.is_model_crs_valid(): + crs = self.data_manager.get_model_crs() + if crs is None or not crs.isValid(): + msg = "Model CRS is not set or invalid. Please select a valid projected CRS in the Model Definition tab." + else: + # Safely get CRS description + try: + crs_desc = crs.description() or crs.authid() or "Unknown" + except Exception: + crs_desc = crs.authid() if hasattr(crs, 'authid') else "Unknown" + msg = f"Model CRS must be projected (in meters), not geographic.\nSelected CRS: {crs_desc}\n\nPlease select a valid projected CRS in the Model Definition tab." + + QMessageBox.critical( + self, + "Invalid Model CRS", + msg, + ) + return # create progress dialog (indeterminate) progress = QProgressDialog("Updating geological model...", "Cancel", 0, 0, self) diff --git a/loopstructural/gui/modelling/model_definition/bounding_box.py b/loopstructural/gui/modelling/model_definition/bounding_box.py index 63a941c..594f70e 100644 --- a/loopstructural/gui/modelling/model_definition/bounding_box.py +++ b/loopstructural/gui/modelling/model_definition/bounding_box.py @@ -2,6 +2,7 @@ import numpy as np from PyQt5.QtWidgets import QWidget +from qgis.core import QgsProject from qgis.PyQt import uic from loopstructural.main.data_manager import default_bounding_box @@ -13,6 +14,8 @@ def __init__(self, parent=None, data_manager=None): super().__init__(parent) ui_path = os.path.join(os.path.dirname(__file__), "bounding_box.ui") uic.loadUi(ui_path, self) + + # Connect bounding box spinbox signals self.originXSpinBox.valueChanged.connect(lambda x: self.onChangeExtent({'xmin': x})) self.maxXSpinBox.valueChanged.connect(lambda x: self.onChangeExtent({'xmax': x})) self.originYSpinBox.valueChanged.connect(lambda y: self.onChangeExtent({'ymin': y})) @@ -21,9 +24,148 @@ def __init__(self, parent=None, data_manager=None): self.maxZSpinBox.valueChanged.connect(lambda z: self.onChangeExtent({'zmax': z})) self.useCurrentViewExtentButton.clicked.connect(self.useCurrentViewExtent) self.selectFromCurrentLayerButton.clicked.connect(self.selectFromCurrentLayer) + + # Connect CRS control signals + self.useProjectCrsRadioButton.toggled.connect(self.onCrsSourceChanged) + self.useCustomCrsRadioButton.toggled.connect(self.onCrsSourceChanged) + self.crsSelector.crsChanged.connect(self.onCrsChanged) + + # Set up callbacks self.data_manager.set_bounding_box_update_callback(self.set_bounding_box) + self.data_manager.set_model_crs_callback(self.update_crs_ui) + + # Initialize CRS UI + self.initialize_crs_ui() self._update_bounding_box_styles() + # Connect to project CRS changes so the widget updates when the project's CRS changes + try: + project = getattr(self.data_manager, 'project', None) or QgsProject.instance() + project.crsChanged.connect(self._onProjectCrsChanged) + except Exception: + # If the signal isn't available or connection fails, ignore to keep widget functional + pass + + def initialize_crs_ui(self): + """Initialize CRS controls with current settings.""" + # Set initial CRS selector value + crs = self.data_manager.get_model_crs() + if crs is not None and crs.isValid(): + self.crsSelector.setCrs(crs) + else: + # Default to project CRS + self.crsSelector.setCrs(self.data_manager.project.crs()) + + # Set radio button based on use_project_crs setting + if self.data_manager._use_project_crs: + self.useProjectCrsRadioButton.setChecked(True) + else: + self.useCustomCrsRadioButton.setChecked(True) + + self.validate_crs() + + def onCrsSourceChanged(self): + """Handle change in CRS source (project vs custom).""" + use_project_crs = self.useProjectCrsRadioButton.isChecked() + self.crsSelector.setEnabled(not use_project_crs) + + if use_project_crs: + # Use project CRS + success, msg = self.data_manager.set_model_crs(None, use_project_crs=True) + else: + # Use custom CRS + crs = self.crsSelector.crs() + success, msg = self.data_manager.set_model_crs(crs, use_project_crs=False) + + self.validate_crs() + + def onCrsChanged(self): + """Handle change in custom CRS selection.""" + if self.useCustomCrsRadioButton.isChecked(): + crs = self.crsSelector.crs() + success, msg = self.data_manager.set_model_crs(crs, use_project_crs=False) + self.validate_crs() + + def update_crs_ui(self, crs, use_project_crs): + """Update UI when model CRS changes externally. + + Parameters + ---------- + crs : QgsCoordinateReferenceSystem or None + The new model CRS + use_project_crs : bool + Whether to use project CRS + """ + # Block signals to avoid recursive updates + self.useProjectCrsRadioButton.blockSignals(True) + self.useCustomCrsRadioButton.blockSignals(True) + self.crsSelector.blockSignals(True) + + try: + if use_project_crs: + self.useProjectCrsRadioButton.setChecked(True) + self.crsSelector.setEnabled(False) + self.crsSelector.setCrs(crs) + + else: + self.useCustomCrsRadioButton.setChecked(True) + self.crsSelector.setEnabled(True) + if crs is not None and crs.isValid(): + self.crsSelector.setCrs(crs) + + self.validate_crs() + finally: + # Unblock signals + self.useProjectCrsRadioButton.blockSignals(False) + self.useCustomCrsRadioButton.blockSignals(False) + self.crsSelector.blockSignals(False) + + def _onProjectCrsChanged(self, crs=None): + """Handle project CRS changes and update UI when the widget is using the project CRS. + + Accept an optional `crs` argument because different QGIS versions may emit the + new CRS or emit no arguments when the project's CRS changes. + """ + # If the signal didn't provide a CRS, try to obtain it from the project's current CRS + if crs is None: + try: + project = getattr(self.data_manager, 'project', None) or QgsProject.instance() + crs = project.crs() + except Exception: + crs = None + + # Only update the UI if the model is configured to use the project CRS + try: + if getattr(self.data_manager, '_use_project_crs', False): + # Update the UI to reflect the new project CRS + self.update_crs_ui(crs, use_project_crs=True) + except Exception: + pass + + def validate_crs(self): + """Validate the selected CRS and update warning label.""" + crs = self.data_manager.get_model_crs() + + if crs is None or not crs.isValid(): + self.crsWarningLabel.setText("⚠ Invalid CRS selected. Model cannot be initialized.") + return False + + if crs.isGeographic(): + # Safely get CRS description + try: + crs_desc = crs.description() or crs.authid() or "Unknown" + except Exception: + crs_desc = crs.authid() if hasattr(crs, 'authid') else "Unknown" + + self.crsWarningLabel.setText( + f"⚠ CRS must be projected (in meters), not geographic.\n" f"Selected: {crs_desc}" + ) + return False + + # CRS is valid and projected + self.crsWarningLabel.setText("") + return True + def set_bounding_box(self, bounding_box): """Populate UI controls with values from a BoundingBox object. diff --git a/loopstructural/gui/modelling/model_definition/bounding_box.ui b/loopstructural/gui/modelling/model_definition/bounding_box.ui index 54e6d89..67523fe 100644 --- a/loopstructural/gui/modelling/model_definition/bounding_box.ui +++ b/loopstructural/gui/modelling/model_definition/bounding_box.ui @@ -162,6 +162,52 @@ + + + + Coordinate Reference System (CRS) + + + + + + Use Project CRS + + + true + + + + + + + Use Custom CRS + + + + + + + false + + + + + + + + + + true + + + color: red; font-weight: bold; + + + + + + @@ -185,6 +231,13 @@ + + + QgsProjectionSelectionWidget + QWidget +
qgis.gui
+ + diff --git a/loopstructural/main/data_manager.py b/loopstructural/main/data_manager.py index 64d816a..fbd4a9b 100644 --- a/loopstructural/main/data_manager.py +++ b/loopstructural/main/data_manager.py @@ -2,7 +2,7 @@ from collections import defaultdict import numpy as np -from qgis.core import QgsPointXY, QgsProject, QgsVectorLayer +from qgis.core import QgsCoordinateReferenceSystem, QgsPointXY, QgsProject, QgsVectorLayer from LoopStructural import FaultTopology, StratigraphicColumn from LoopStructural.datatypes import BoundingBox @@ -69,6 +69,9 @@ def __init__(self, *, project=None, mapCanvas=None, logger=None): self.dem_callback = None self.widget_settings = {} self.feature_data = defaultdict(dict) + self._model_crs = None + self._use_project_crs = True + self.model_crs_callback = None def onSaveProject(self): """Save project data.""" @@ -171,6 +174,100 @@ def get_bounding_box(self): """Get the current bounding box.""" return self._bounding_box + def set_model_crs(self, crs, use_project_crs=False): + """Set the model CRS. + + Parameters + ---------- + crs : QgsCoordinateReferenceSystem or None + The CRS to use for the model. If None and use_project_crs is True, + will use the project CRS. + use_project_crs : bool + If True, use the project CRS instead of a custom CRS. + + Returns + ------- + tuple + (success: bool, message: str) + """ + self._use_project_crs = use_project_crs + + if use_project_crs: + crs = self.project.crs() + + # Validate CRS + if crs is None or not crs.isValid(): + self._model_crs = None + msg = "Model CRS is not valid." + self.logger(message=msg, log_level=2) + if self.model_crs_callback: + self.model_crs_callback(self._model_crs, self._use_project_crs) + return False, msg + + # Check if CRS is projected (not geographic) + if crs.isGeographic(): + self._model_crs = None + # Safely get CRS description + try: + crs_desc = crs.description() or crs.authid() or "Unknown" + except Exception: + crs_desc = crs.authid() if hasattr(crs, 'authid') else "Unknown" + msg = f"Model CRS must be projected (in meters), not geographic. Selected CRS: {crs_desc}" + self.logger(message=msg, log_level=2) + if self.model_crs_callback: + self.model_crs_callback(self._model_crs, self._use_project_crs) + return False, msg + + self._model_crs = crs + # Safely get CRS description + try: + crs_desc = crs.description() or "Unknown" + crs_id = crs.authid() or "Unknown" + except Exception: + crs_desc = "Unknown" + crs_id = crs.authid() if hasattr(crs, 'authid') else "Unknown" + msg = f"Model CRS set to: {crs_desc} ({crs_id})" + self.logger(message=msg, log_level=3) + + if self.model_crs_callback: + self.model_crs_callback(self._model_crs, self._use_project_crs) + + return True, msg + + def get_model_crs(self): + """Get the model CRS. + + Returns + ------- + QgsCoordinateReferenceSystem or None + The model CRS, or None if not set. + """ + if self._use_project_crs: + return self.project.crs() + return self._model_crs + + def is_model_crs_valid(self): + """Check if the model CRS is valid and projected. + + Returns + ------- + bool + True if the model CRS is valid and projected, False otherwise. + """ + crs = self.get_model_crs() + if crs is None or not crs.isValid(): + return False + if crs.isGeographic(): + return False + return True + + def set_model_crs_callback(self, callback): + """Set the callback for when the model CRS is updated.""" + self.model_crs_callback = callback + # Trigger callback with current values + if self.model_crs_callback: + self.model_crs_callback(self.get_model_crs(), self._use_project_crs) + def set_elevation(self, elevation): """Set the elevation for the model.""" self.elevation = elevation @@ -378,15 +475,16 @@ def update_stratigraphy(self): """Update the foliation features in the model manager.""" print("Updating stratigraphy...") if self._model_manager is not None: + model_crs = self.get_model_crs() if self._basal_contacts is not None: self._model_manager.update_contact_traces( - qgsLayerToGeoDataFrame(self._basal_contacts['layer']), + qgsLayerToGeoDataFrame(self._basal_contacts['layer'], target_crs=model_crs), unit_name_field=self._basal_contacts['unitname_field'], ) if self._structural_orientations is not None: print("Updating structural orientations...") self._model_manager.update_structural_data( - qgsLayerToGeoDataFrame(self._structural_orientations['layer']), + qgsLayerToGeoDataFrame(self._structural_orientations['layer'], target_crs=model_crs), strike_field=self._structural_orientations['strike_field'], dip_field=self._structural_orientations['dip_field'], unit_name_field=self._structural_orientations['unitname_field'], @@ -412,8 +510,9 @@ def update_faults(self): self._fault_topology.remove_fault(fault) self.fault_adjacency = np.zeros((len(unique_faults), len(unique_faults)), dtype=int) if self._model_manager is not None: + model_crs = self.get_model_crs() self._model_manager.update_fault_points( - qgsLayerToGeoDataFrame(self._fault_traces['layer']), + qgsLayerToGeoDataFrame(self._fault_traces['layer'], target_crs=model_crs), fault_name_field=self._fault_traces['fault_name_field'], fault_dip_field=self._fault_traces['fault_dip_field'], fault_pitch_field=self._fault_traces.get('fault_pitch_field', None), @@ -437,6 +536,22 @@ def clear_data(self): self._fault_traces = None self._structural_orientations = None + def _get_model_crs_authid(self): + """Get the model CRS authid string for serialization. + + Returns + ------- + str or None + CRS authid string (e.g., 'EPSG:32633') or None if CRS is not valid. + """ + if not self._model_crs: + return None + if not isinstance(self._model_crs, QgsCoordinateReferenceSystem): + return None + if not self._model_crs.isValid(): + return None + return self._model_crs.authid() + def to_dict(self): """Convert the data manager to a dictionary.""" # Create copies of the dictionaries to avoid modifying the originals @@ -476,6 +591,8 @@ def to_dict(self): 'use_dem': self.use_dem, 'elevation': self.elevation, 'widget_settings': self.widget_settings, + 'model_crs': self._get_model_crs_authid(), + 'use_project_crs': self._use_project_crs, } def from_dict(self, data): @@ -514,6 +631,14 @@ def from_dict(self, data): self.stratigraphic_column_callback() if 'widget_settings' in data: self.widget_settings = data['widget_settings'] + + # Load model CRS settings + if 'use_project_crs' in data: + self._use_project_crs = data['use_project_crs'] + if 'model_crs' in data and data['model_crs'] is not None: + crs = QgsCoordinateReferenceSystem(data['model_crs']) + if crs.isValid(): + self.set_model_crs(crs, use_project_crs=self._use_project_crs) def update_from_dict(self, data): """Update the data manager from a dictionary.""" @@ -590,6 +715,17 @@ def update_from_dict(self, data): self.widget_settings = data['widget_settings'] else: self.widget_settings = {} + + # Load model CRS settings + if 'use_project_crs' in data: + self._use_project_crs = data['use_project_crs'] + else: + self._use_project_crs = True + + if 'model_crs' in data and data['model_crs'] is not None: + crs = QgsCoordinateReferenceSystem(data['model_crs']) + if crs.isValid(): + self.set_model_crs(crs, use_project_crs=self._use_project_crs) if self.stratigraphic_column_callback: self.stratigraphic_column_callback() @@ -644,9 +780,10 @@ def add_foliation_to_model(self, foliation_name: str, *, folded_feature_name=Non if foliation_name not in self.feature_data: raise ValueError(f"Foliation '{foliation_name}' does not exist in the data manager.") foliation_data = self.feature_data[foliation_name] + model_crs = self.get_model_crs() for layer in foliation_data.values(): layer['df'] = qgsLayerToGeoDataFrame( - layer['layer'] + layer['layer'], target_crs=model_crs ) # Convert QgsVectorLayer to GeoDataFrame if self._model_manager: self._model_manager.add_foliation( diff --git a/loopstructural/main/vectorLayerWrapper.py b/loopstructural/main/vectorLayerWrapper.py index b4ed3e9..f744d0c 100644 --- a/loopstructural/main/vectorLayerWrapper.py +++ b/loopstructural/main/vectorLayerWrapper.py @@ -1,5 +1,7 @@ +import logging import os import tempfile +from typing import Optional import geopandas as gpd import numpy as np @@ -27,6 +29,8 @@ from qgis.PyQt.QtCore import QDateTime, QVariant from shapely.geometry import LineString, MultiLineString, MultiPoint, MultiPolygon, Point, Polygon +logger = logging.getLogger(__name__) + def qgsRasterToGdalDataset(rlayer: QgsRasterLayer): """ @@ -88,28 +92,102 @@ def qgsRasterToGdalDataset(rlayer: QgsRasterLayer): return ds -def qgsLayerToGeoDataFrame(layer) -> gpd.GeoDataFrame: +def _get_crs_id(crs): + """Get a safe string identifier for a CRS. + + Parameters + ---------- + crs : QgsCoordinateReferenceSystem or None + The CRS to get an ID for + + Returns + ------- + str + CRS authid or "Unknown" if unavailable + """ + if crs and crs.isValid(): + try: + return crs.authid() or "Unknown" + except Exception: + return "Unknown" + return "Unknown" + + +def qgsLayerToGeoDataFrame(layer, target_crs=None) -> Optional[gpd.GeoDataFrame]: + """Convert a QgsVectorLayer to a GeoDataFrame, optionally transforming to a target CRS. + + Parameters + ---------- + layer : QgsVectorLayer + The vector layer to convert + target_crs : QgsCoordinateReferenceSystem, optional + If provided, all geometries will be transformed to this CRS. + If None, the layer's source CRS is used. + + Returns + ------- + gpd.GeoDataFrame + GeoDataFrame with geometries in the specified CRS + """ if layer is None: return None + features = layer.getFeatures() fields = layer.fields() data = {'geometry': []} for f in fields: data[f.name()] = [] + + # Set up coordinate transformation if needed + transform = None + source_crs = layer.sourceCrs() + output_crs = source_crs + + if target_crs is not None and target_crs.isValid(): + if source_crs.isValid() and source_crs != target_crs: + transform = QgsCoordinateTransform(source_crs, target_crs, QgsProject.instance()) + output_crs = target_crs + for feature in features: geom = feature.geometry() if geom.isEmpty(): continue - data['geometry'].append(geom) + + # Transform geometry if needed + if transform is not None: + geom_copy = QgsGeometry(geom) + try: + result = geom_copy.transform(transform) + if result != 0: + # Transform returned error code + logger.warning( + f"Failed to transform geometry (error code {result}). " + f"Source CRS: {_get_crs_id(source_crs)}, Target CRS: {_get_crs_id(target_crs)}. " + f"Skipping feature." + ) + continue + except Exception as e: + # If transformation fails, log warning and skip this feature + logger.exception( + f"Exception during CRS transformation: {e}. " + f"Source CRS: {_get_crs_id(source_crs)}, Target CRS: {_get_crs_id(target_crs)}. " + f"Skipping feature." + ) + continue + else: + data['geometry'].append(geom) + + # Copy field values 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()) + + return gpd.GeoDataFrame(data, crs=output_crs.authid()) -def qgsLayerToDataFrame(src, dtm=None) -> pd.DataFrame: +def qgsLayerToDataFrame(src, dtm=None) -> Optional[pd.DataFrame]: """ Convert a vector layer or processing feature source to a pandas DataFrame. Samples geometry using points or vertices of lines/polygons. diff --git a/tests/qgis/test_model_crs.py b/tests/qgis/test_model_crs.py new file mode 100644 index 0000000..a099709 --- /dev/null +++ b/tests/qgis/test_model_crs.py @@ -0,0 +1,127 @@ +import unittest +from unittest.mock import MagicMock, Mock + +from qgis.core import QgsCoordinateReferenceSystem, QgsProject + +from loopstructural.main.data_manager import ModellingDataManager + + +class TestModelCRS(unittest.TestCase): + """Unit tests for Model CRS functionality.""" + + def setUp(self): + """Set up test fixtures.""" + # Create mock objects + self.mock_project = Mock(spec=QgsProject) + self.mock_canvas = Mock() + self.mock_logger = Mock() + + # Set up mock project CRS (WGS84 UTM 55S - projected) + self.mock_project_crs = QgsCoordinateReferenceSystem("EPSG:32755") + self.mock_project.crs.return_value = self.mock_project_crs + + # Initialize data manager + self.data_manager = ModellingDataManager( + project=self.mock_project, + mapCanvas=self.mock_canvas, + logger=self.mock_logger + ) + + def test_initial_crs_uses_project_crs_flag(self): + """Test that initial state has use_project_crs=True.""" + self.assertTrue(self.data_manager._use_project_crs) + + def test_set_valid_projected_crs(self): + """Test setting a valid projected CRS.""" + # Create a projected CRS (WGS84 / UTM zone 33N) + crs = QgsCoordinateReferenceSystem("EPSG:32633") + + success, msg = self.data_manager.set_model_crs(crs, use_project_crs=False) + + self.assertTrue(success) + self.assertEqual(self.data_manager.get_model_crs(), crs) + self.assertFalse(self.data_manager._use_project_crs) + + def test_reject_geographic_crs(self): + """Test that geographic CRS is rejected.""" + # Create a geographic CRS (WGS84) + crs = QgsCoordinateReferenceSystem("EPSG:4326") + + success, msg = self.data_manager.set_model_crs(crs, use_project_crs=False) + + self.assertFalse(success) + self.assertIn("geographic", msg.lower()) + self.assertIsNone(self.data_manager._model_crs) + + def test_reject_invalid_crs(self): + """Test that invalid CRS is rejected.""" + # Create an invalid CRS + crs = QgsCoordinateReferenceSystem() + + success, msg = self.data_manager.set_model_crs(crs, use_project_crs=False) + + self.assertFalse(success) + self.assertIn("not valid", msg.lower()) + self.assertIsNone(self.data_manager._model_crs) + + def test_use_project_crs(self): + """Test using project CRS.""" + success, msg = self.data_manager.set_model_crs(None, use_project_crs=True) + + self.assertTrue(success) + self.assertTrue(self.data_manager._use_project_crs) + # Should return project CRS when use_project_crs is True + self.assertEqual(self.data_manager.get_model_crs(), self.mock_project_crs) + + def test_is_model_crs_valid_with_projected(self): + """Test CRS validation with a projected CRS.""" + crs = QgsCoordinateReferenceSystem("EPSG:32633") + self.data_manager.set_model_crs(crs, use_project_crs=False) + + self.assertTrue(self.data_manager.is_model_crs_valid()) + + def test_is_model_crs_valid_with_geographic(self): + """Test CRS validation rejects geographic CRS.""" + crs = QgsCoordinateReferenceSystem("EPSG:4326") + self.data_manager.set_model_crs(crs, use_project_crs=False) + + self.assertFalse(self.data_manager.is_model_crs_valid()) + + def test_crs_persistence_to_dict(self): + """Test CRS is saved in to_dict.""" + crs = QgsCoordinateReferenceSystem("EPSG:32633") + self.data_manager.set_model_crs(crs, use_project_crs=False) + + data_dict = self.data_manager.to_dict() + + self.assertEqual(data_dict['model_crs'], 'EPSG:32633') + self.assertFalse(data_dict['use_project_crs']) + + def test_crs_persistence_from_dict(self): + """Test CRS is restored from from_dict.""" + # Create a data dict with CRS info + data = { + 'model_crs': 'EPSG:32633', + 'use_project_crs': False, + } + + # Clear current CRS + self.data_manager._model_crs = None + self.data_manager._use_project_crs = True + + # Load from dict (just the CRS part - simplified test) + if 'use_project_crs' in data: + self.data_manager._use_project_crs = data['use_project_crs'] + if 'model_crs' in data and data['model_crs'] is not None: + crs = QgsCoordinateReferenceSystem(data['model_crs']) + if crs.isValid(): + self.data_manager.set_model_crs(crs, use_project_crs=self.data_manager._use_project_crs) + + # Verify + self.assertFalse(self.data_manager._use_project_crs) + restored_crs = self.data_manager.get_model_crs() + self.assertEqual(restored_crs.authid(), 'EPSG:32633') + + +if __name__ == '__main__': + unittest.main()