profile_interpreter.py (8971B)
1 """Profile Interpreter — proof of concept. 2 3 Captures mouse clicks on the QGIS Elevation Profile canvas and places PointZ 4 features at the corresponding along-section position and elevation. This lets 5 you digitize interpretation points (geological boundaries, borehole picks, 6 geophysical layer tops, ...) directly on a profile view. 7 8 The conversion path is: 9 10 click position (canvas pixels) 11 -> QgsElevationProfileCanvas.canvasPointToPlotPoint() 12 -> QgsProfilePoint(distance_along_section, elevation) 13 -> interpolate XY along the profile curve at `distance` 14 -> QgsPoint(x, y, z = elevation) -> PointZ feature 15 """ 16 17 import math 18 import os 19 20 from qgis.PyQt.QtCore import Qt, QPointF 21 from qgis.PyQt.QtGui import QIcon 22 from qgis.PyQt.QtWidgets import QAction 23 24 try: 25 _LEFT_BUTTON = Qt.MouseButton.LeftButton # PyQt6 / QGIS 4 26 except AttributeError: 27 _LEFT_BUTTON = Qt.LeftButton # PyQt5 / QGIS 3 28 29 from qgis.core import ( 30 Qgis, 31 QgsCoordinateTransform, 32 QgsFeature, 33 QgsField, 34 QgsGeometry, 35 QgsMessageLog, 36 QgsPoint, 37 QgsPointXY, 38 QgsProject, 39 QgsVectorDataProvider, 40 QgsVectorLayer, 41 QgsWkbTypes, 42 ) 43 from qgis.PyQt.QtCore import QMetaType 44 from qgis.gui import QgsElevationProfileCanvas, QgsPlotTool 45 46 try: 47 from qgis.gui import QgsPlotToolPan 48 except ImportError: # older builds expose it under a different path 49 QgsPlotToolPan = None 50 51 try: 52 _POINT_GEOMETRY = QgsWkbTypes.GeometryType.PointGeometry # PyQt6 / QGIS 4 53 _ADD_FEATURES_CAP = QgsVectorDataProvider.Capability.AddFeatures 54 except AttributeError: 55 _POINT_GEOMETRY = QgsWkbTypes.PointGeometry # PyQt5 / QGIS 3 56 _ADD_FEATURES_CAP = QgsVectorDataProvider.AddFeatures 57 58 59 MENU = 'Profile Interpreter' 60 LAYER_NAME = 'Profile interpretations' 61 62 63 def _is_suitable_target(layer): 64 """True if layer is a writable PointZ vector layer.""" 65 if layer is None or not isinstance(layer, QgsVectorLayer): 66 return False 67 if layer.geometryType() != _POINT_GEOMETRY: 68 return False 69 if not QgsWkbTypes.hasZ(layer.wkbType()): 70 return False 71 caps = layer.dataProvider().capabilities() 72 return bool(caps & _ADD_FEATURES_CAP) 73 74 75 def _point_attributes(field_names, feature_id, distance, elevation): 76 """Return {field: value} for only the names present in field_names.""" 77 values = {'id': feature_id, 'distance': distance, 'elevation': elevation} 78 return {k: v for k, v in values.items() if k in field_names} 79 80 81 def _project_xy(xy, src_crs, dst_crs): 82 """Map x/y from src_crs to dst_crs; identity if either CRS is invalid or they are equal.""" 83 if (not src_crs.isValid() or not dst_crs.isValid() or 84 src_crs == dst_crs): 85 return xy.x(), xy.y() 86 ct = QgsCoordinateTransform(src_crs, dst_crs, QgsProject.instance()) 87 p = ct.transform(QgsPointXY(xy.x(), xy.y())) 88 return p.x(), p.y() 89 90 91 def _next_feature_id(layer): 92 """Next integer id for layer: max existing 'id' + 1, or 1. 93 Returns None if the layer has no 'id' field.""" 94 names = {f.name() for f in layer.fields()} 95 if 'id' not in names: 96 return None 97 max_id = 0 98 for feat in layer.getFeatures(): 99 try: 100 val = feat['id'] 101 except (KeyError, IndexError): 102 continue 103 if val is not None and int(val) > max_id: 104 max_id = int(val) 105 return max_id + 1 106 107 108 class _ProfilePickTool(QgsPlotTool): 109 """A plot tool that forwards left-button releases to a callback.""" 110 111 def __init__(self, canvas, on_pick): 112 super().__init__(canvas, 'Interpret on Profile') 113 self._on_pick = on_pick 114 115 def plotReleaseEvent(self, event): 116 if event.button() != _LEFT_BUTTON: 117 super().plotReleaseEvent(event) 118 return 119 self._on_pick(event) 120 121 122 class ProfileInterpreterPlugin: 123 def __init__(self, iface): 124 self.iface = iface 125 self._action = None 126 self._canvas = None 127 self._pick_tool = None 128 self._pan_tool = None 129 self._layer = None 130 131 def initGui(self): 132 icon_path = os.path.join(os.path.dirname(__file__), 'icon.svg') 133 self._action = QAction(QIcon(icon_path), 'Interpret on Profile', self.iface.mainWindow()) 134 self._action.setCheckable(True) 135 self._action.toggled.connect(self._set_active) 136 self.iface.addPluginToMenu(MENU, self._action) 137 self.iface.addToolBarIcon(self._action) 138 139 def unload(self): 140 self._set_active(False) 141 self.iface.removeToolBarIcon(self._action) 142 self.iface.removePluginMenu(MENU, self._action) 143 self._action = None 144 145 def _set_active(self, active): 146 if not active: 147 self._deactivate() 148 return 149 150 canvas = self._find_profile_canvas() 151 if canvas is None: 152 self.iface.messageBar().pushMessage( 153 MENU, 154 'Open the Elevation Profile panel first ' 155 '(View > Elevation Profile), then enable this tool.', 156 level=Qgis.Warning, 157 ) 158 if self._action is not None: 159 self._action.setChecked(False) 160 return 161 162 self._canvas = canvas 163 self._pick_tool = _ProfilePickTool(canvas, self._on_pick) 164 canvas.setTool(self._pick_tool) 165 self.iface.messageBar().pushMessage( 166 MENU, 167 'Click on the elevation profile to place interpretation points.', 168 level=Qgis.Info, 169 ) 170 171 def _deactivate(self): 172 if self._canvas is not None and self._pick_tool is not None: 173 if QgsPlotToolPan is not None: 174 self._pan_tool = QgsPlotToolPan(self._canvas) 175 else: 176 self._pan_tool = QgsPlotTool(self._canvas, 'Inactive') 177 self._canvas.setTool(self._pan_tool) 178 self._pick_tool = None 179 self._canvas = None 180 181 def _find_profile_canvas(self): 182 main_window = self.iface.mainWindow() 183 canvases = main_window.findChildren(QgsElevationProfileCanvas) 184 for canvas in canvases: 185 if canvas.isVisible(): 186 return canvas 187 return canvases[0] if canvases else None 188 189 def _target_layer(self, crs): 190 active = self.iface.activeLayer() 191 if _is_suitable_target(active): 192 return active 193 return self._interpretation_layer(crs) 194 195 def _on_pick(self, event): 196 canvas = self._canvas 197 profile_point = canvas.canvasPointToPlotPoint(QPointF(event.pos())) 198 distance = profile_point.distance() 199 elevation = profile_point.elevation() 200 if math.isnan(distance) or math.isnan(elevation): 201 return 202 203 curve = canvas.profileCurve() 204 if curve is None: 205 return 206 point_geom = QgsGeometry(curve.clone()).interpolate(distance) 207 if point_geom.isEmpty(): 208 return 209 xy = point_geom.asPoint() 210 211 layer = self._target_layer(canvas.crs()) 212 x, y = _project_xy(xy, canvas.crs(), layer.crs()) 213 feature = QgsFeature(layer.fields()) 214 feature.setGeometry(QgsGeometry(QgsPoint(x, y, elevation))) 215 216 feature_id = _next_feature_id(layer) 217 attrs = _point_attributes( 218 {f.name() for f in layer.fields()}, 219 feature_id, distance, elevation, 220 ) 221 for k, v in attrs.items(): 222 feature[k] = v 223 224 ok, _ = layer.dataProvider().addFeatures([feature]) 225 if not ok: 226 self.iface.messageBar().pushMessage( 227 MENU, 228 'Failed to add feature to layer.', 229 level=Qgis.Warning, 230 ) 231 return 232 layer.updateExtents() 233 layer.triggerRepaint() 234 235 if hasattr(canvas, 'refresh'): 236 canvas.refresh() 237 238 id_label = f' #{feature_id}' if feature_id is not None else '' 239 QgsMessageLog.logMessage( 240 f'placed point{id_label} at x={x:.3f} y={y:.3f} z={elevation:.3f} ' 241 f'(distance={distance:.3f})', 242 MENU, level=Qgis.Info, 243 ) 244 self.iface.messageBar().pushMessage( 245 MENU, 246 f'Point at x={xy.x():.1f}, y={xy.y():.1f}, ' 247 f'elevation={elevation:.2f}', 248 level=Qgis.Success, 249 duration=2, 250 ) 251 252 def _interpretation_layer(self, crs): 253 project = QgsProject.instance() 254 if self._layer is not None and project.mapLayer(self._layer.id()): 255 return self._layer 256 257 uri = f'PointZ?crs={crs.authid()}' if crs.isValid() else 'PointZ' 258 layer = QgsVectorLayer(uri, LAYER_NAME, 'memory') 259 provider = layer.dataProvider() 260 provider.addAttributes([ 261 QgsField('id', QMetaType.Type.Int), 262 QgsField('distance', QMetaType.Type.Double), 263 QgsField('elevation', QMetaType.Type.Double), 264 ]) 265 layer.updateFields() 266 project.addMapLayer(layer) 267 self._layer = layer 268 return layer