profile_interpreter.py (22903B)
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, QToolBar 23 24 try: 25 _LEFT_BUTTON = Qt.MouseButton.LeftButton # PyQt6 / QGIS 4 26 _RIGHT_BUTTON = Qt.MouseButton.RightButton 27 except AttributeError: 28 _LEFT_BUTTON = Qt.LeftButton # PyQt5 / QGIS 3 29 _RIGHT_BUTTON = Qt.RightButton 30 31 from qgis.core import ( 32 Qgis, 33 QgsApplication, 34 QgsCoordinateTransform, 35 QgsFeature, 36 QgsField, 37 QgsGeometry, 38 QgsMessageLog, 39 QgsPoint, 40 QgsPointXY, 41 QgsProfilePoint, 42 QgsProject, 43 QgsVectorDataProvider, 44 QgsVectorLayer, 45 QgsWkbTypes, 46 ) 47 from qgis.PyQt.QtCore import QMetaType 48 from qgis.gui import QgsElevationProfileCanvas, QgsPlotTool 49 50 try: 51 from qgis.gui import QgsPlotToolPan 52 except ImportError: # older builds expose it under a different path 53 QgsPlotToolPan = None 54 55 try: 56 _POINT_GEOMETRY = QgsWkbTypes.GeometryType.PointGeometry # PyQt6 / QGIS 4 57 _ADD_FEATURES_CAP = QgsVectorDataProvider.Capability.AddFeatures 58 _DELETE_FEATURES_CAP = QgsVectorDataProvider.Capability.DeleteFeatures 59 _CHANGE_GEOMETRIES_CAP = QgsVectorDataProvider.Capability.ChangeGeometries 60 _CHANGE_ATTRIBUTES_CAP = QgsVectorDataProvider.Capability.ChangeAttributeValues 61 except AttributeError: 62 _POINT_GEOMETRY = QgsWkbTypes.PointGeometry # PyQt5 / QGIS 3 63 _ADD_FEATURES_CAP = QgsVectorDataProvider.AddFeatures 64 _DELETE_FEATURES_CAP = QgsVectorDataProvider.DeleteFeatures 65 _CHANGE_GEOMETRIES_CAP = QgsVectorDataProvider.ChangeGeometries 66 _CHANGE_ATTRIBUTES_CAP = QgsVectorDataProvider.ChangeAttributeValues 67 68 69 MENU = 'Profile Interpreter' 70 LAYER_NAME = 'Profile interpretations' 71 HIT_TOLERANCE_PX = 10.0 72 73 74 def _is_suitable_target(layer): 75 """True if layer is a writable PointZ vector layer.""" 76 if layer is None or not isinstance(layer, QgsVectorLayer): 77 return False 78 if layer.geometryType() != _POINT_GEOMETRY: 79 return False 80 if not QgsWkbTypes.hasZ(layer.wkbType()): 81 return False 82 caps = layer.dataProvider().capabilities() 83 return bool(caps & _ADD_FEATURES_CAP) 84 85 86 def _point_attributes(field_names, feature_id, distance, elevation): 87 """Return {field: value} for only the names present in field_names.""" 88 values = {'id': feature_id, 'distance': distance, 'elevation': elevation} 89 return {k: v for k, v in values.items() if k in field_names} 90 91 92 def _project_xy(xy, src_crs, dst_crs): 93 """Map x/y from src_crs to dst_crs; identity if either CRS is invalid or they are equal.""" 94 if (not src_crs.isValid() or not dst_crs.isValid() or 95 src_crs == dst_crs): 96 return xy.x(), xy.y() 97 ct = QgsCoordinateTransform(src_crs, dst_crs, QgsProject.instance()) 98 p = ct.transform(QgsPointXY(xy.x(), xy.y())) 99 return p.x(), p.y() 100 101 102 def _next_feature_id(layer): 103 """Next integer id for layer: max existing 'id' + 1, or 1. 104 Returns None if the layer has no 'id' field.""" 105 names = {f.name() for f in layer.fields()} 106 if 'id' not in names: 107 return None 108 max_id = 0 109 for feat in layer.getFeatures(): 110 try: 111 val = feat['id'] 112 except (KeyError, IndexError): 113 continue 114 if val is not None and int(val) > max_id: 115 max_id = int(val) 116 return max_id + 1 117 118 119 def _can_delete(layer): 120 """True if features can be removed from layer's provider.""" 121 if layer is None: 122 return False 123 caps = layer.dataProvider().capabilities() 124 return bool(caps & _DELETE_FEATURES_CAP) 125 126 127 def _can_move(layer): 128 """True if feature geometries can be changed in layer's provider.""" 129 if layer is None: 130 return False 131 caps = layer.dataProvider().capabilities() 132 return bool(caps & _CHANGE_GEOMETRIES_CAP) 133 134 135 def _nearest_hit(candidates, click_x, click_y, tolerance_px): 136 """Feature id of the (fid, px, py) candidate nearest the click within tolerance, or None.""" 137 best_fid = None 138 best_dist = tolerance_px 139 for fid, px, py in candidates: 140 dist = math.hypot(px - click_x, py - click_y) 141 if dist <= best_dist: 142 best_fid = fid 143 best_dist = dist 144 return best_fid 145 146 147 def _added_feature_id(added, fallback): 148 """Provider-assigned id of the feature just written, else the fallback feature's id.""" 149 feature = added[0] if added else fallback 150 return feature.id() 151 152 153 def _find_panel_toolbar(canvas, main_window): 154 """First QToolBar in the panel holding canvas, or None. 155 156 Walks up from the canvas and stops at the main window so the plugin never 157 lands its action on one of the main QGIS toolbars. 158 """ 159 widget = canvas.parentWidget() 160 while widget is not None and widget is not main_window: 161 bars = widget.findChildren(QToolBar) 162 if bars: 163 return bars[0] 164 widget = widget.parentWidget() 165 return None 166 167 168 class _ProfilePickTool(QgsPlotTool): 169 """A plot tool that forwards left-button releases to a callback. 170 171 Right-button releases call on_cancel, if given, before default handling. 172 """ 173 174 def __init__(self, canvas, on_pick, on_cancel=None): 175 super().__init__(canvas, 'Interpret on Profile') 176 self._on_pick = on_pick 177 self._on_cancel = on_cancel 178 179 def plotReleaseEvent(self, event): 180 if event.button() == _LEFT_BUTTON: 181 self._on_pick(event) 182 return 183 if event.button() == _RIGHT_BUTTON and self._on_cancel is not None: 184 self._on_cancel() 185 super().plotReleaseEvent(event) 186 187 188 class ProfileInterpreterPlugin: 189 def __init__(self, iface): 190 self.iface = iface 191 self._action = None 192 self._canvas = None 193 self._pick_tool = None 194 self._pan_tool = None 195 self._layer = None 196 self._undo_action = None 197 self._edit_action = None 198 self._delete_action = None 199 self._panel_toolbar = None 200 self._history = [] 201 self._mode = 'add' 202 self._selected = None 203 204 def initGui(self): 205 icon_path = os.path.join(os.path.dirname(__file__), 'icon.svg') 206 self._action = QAction(QIcon(icon_path), 'Interpret on Profile', self.iface.mainWindow()) 207 self._action.setCheckable(True) 208 self._action.toggled.connect(self._set_active) 209 self.iface.addPluginToMenu(MENU, self._action) 210 self.iface.addToolBarIcon(self._action) 211 212 self._undo_action = QAction( 213 QgsApplication.getThemeIcon('mActionUndo.svg'), 214 'Undo Point', self.iface.mainWindow()) 215 self._undo_action.setToolTip('Remove the last interpretation point') 216 self._undo_action.setEnabled(False) 217 self._undo_action.triggered.connect(self._undo_last_point) 218 self.iface.addPluginToMenu(MENU, self._undo_action) 219 220 self._edit_action = QAction( 221 QgsApplication.getThemeIcon('mActionVertexTool.svg'), 222 'Edit Point', self.iface.mainWindow()) 223 self._edit_action.setToolTip( 224 'Click a point to select it, then click its new position') 225 self._edit_action.setCheckable(True) 226 self._edit_action.setEnabled(False) 227 self._edit_action.toggled.connect(lambda on: self._set_mode('edit', on)) 228 self.iface.addPluginToMenu(MENU, self._edit_action) 229 230 self._delete_action = QAction( 231 QgsApplication.getThemeIcon('mActionDeleteSelectedFeatures.svg'), 232 'Delete Point', self.iface.mainWindow()) 233 self._delete_action.setToolTip('Click a point to delete it') 234 self._delete_action.setCheckable(True) 235 self._delete_action.setEnabled(False) 236 self._delete_action.toggled.connect(lambda on: self._set_mode('delete', on)) 237 self.iface.addPluginToMenu(MENU, self._delete_action) 238 239 def unload(self): 240 self._set_active(False) 241 self.iface.removeToolBarIcon(self._action) 242 for action in (self._action, self._undo_action, self._edit_action, 243 self._delete_action): 244 self.iface.removePluginMenu(MENU, action) 245 self._action = None 246 self._undo_action = None 247 self._edit_action = None 248 self._delete_action = None 249 self._history.clear() 250 251 def _set_active(self, active): 252 if not active: 253 self._deactivate() 254 return 255 256 canvas = self._find_profile_canvas() 257 if canvas is None: 258 self.iface.messageBar().pushMessage( 259 MENU, 260 'Open the Elevation Profile panel first ' 261 '(View > Elevation Profile), then enable this tool.', 262 level=Qgis.Warning, 263 ) 264 if self._action is not None: 265 self._action.setChecked(False) 266 return 267 268 self._canvas = canvas 269 self._pick_tool = _ProfilePickTool(canvas, self._on_pick, self._clear_selection) 270 canvas.setTool(self._pick_tool) 271 self._attach_panel_actions(canvas) 272 self.iface.messageBar().pushMessage( 273 MENU, 274 'Click on the elevation profile to place interpretation points.', 275 level=Qgis.Info, 276 ) 277 278 def _deactivate(self): 279 self._detach_panel_actions() 280 if self._canvas is not None and self._pick_tool is not None: 281 if QgsPlotToolPan is not None: 282 self._pan_tool = QgsPlotToolPan(self._canvas) 283 else: 284 self._pan_tool = QgsPlotTool(self._canvas, 'Inactive') 285 self._canvas.setTool(self._pan_tool) 286 self._pick_tool = None 287 self._canvas = None 288 289 def _panel_actions(self): 290 actions = (self._undo_action, self._edit_action, self._delete_action) 291 return [a for a in actions if a is not None] 292 293 def _attach_panel_actions(self, canvas): 294 """Enable the mode actions and put all panel actions on the profile panel toolbar.""" 295 for action in (self._edit_action, self._delete_action): 296 if action is not None: 297 action.setEnabled(True) 298 actions = self._panel_actions() 299 if not actions: 300 return 301 toolbar = _find_panel_toolbar(canvas, self.iface.mainWindow()) 302 if toolbar is None: 303 return 304 for action in actions: 305 toolbar.addAction(action) 306 self._panel_toolbar = toolbar 307 308 def _detach_panel_actions(self): 309 for action in (self._edit_action, self._delete_action): 310 if action is not None: 311 action.setChecked(False) 312 action.setEnabled(False) 313 self._clear_selection() 314 self._mode = 'add' 315 if self._panel_toolbar is not None: 316 for action in self._panel_actions(): 317 self._panel_toolbar.removeAction(action) 318 self._panel_toolbar = None 319 320 def _set_mode(self, mode, checked): 321 """Switch between add, edit and delete; edit and delete are mutually exclusive.""" 322 self._clear_selection() 323 if not checked: 324 if self._mode == mode: 325 self._mode = 'add' 326 return 327 self._mode = mode 328 other = self._delete_action if mode == 'edit' else self._edit_action 329 if other is not None: 330 other.setChecked(False) 331 hint = ('Click a point to select it, then click its new position.' 332 if mode == 'edit' else 'Click a point to delete it.') 333 self.iface.messageBar().pushMessage(MENU, hint, level=Qgis.Info, duration=3) 334 335 def _clear_selection(self): 336 if self._selected is not None: 337 layer, _ = self._selected 338 if QgsProject.instance().mapLayer(layer.id()) is not None: 339 layer.removeSelection() 340 self._selected = None 341 342 def _refresh_undo_state(self): 343 if self._undo_action is not None: 344 self._undo_action.setEnabled(bool(self._history)) 345 346 def _refresh_layer(self, layer): 347 layer.updateExtents() 348 layer.triggerRepaint() 349 if self._canvas is not None and hasattr(self._canvas, 'refresh'): 350 self._canvas.refresh() 351 352 def _undo_last_point(self): 353 """Delete the most recently placed point, skipping layers since removed.""" 354 project = QgsProject.instance() 355 while self._history: 356 layer, feature_id = self._history.pop() 357 if project.mapLayer(layer.id()) is None: 358 continue 359 if not _can_delete(layer): 360 self.iface.messageBar().pushMessage( 361 MENU, 362 f'Layer "{layer.name()}" does not support deleting features.', 363 level=Qgis.Warning, 364 ) 365 break 366 if not layer.dataProvider().deleteFeatures([feature_id]): 367 self.iface.messageBar().pushMessage( 368 MENU, 369 'Failed to remove the last point.', 370 level=Qgis.Warning, 371 ) 372 break 373 self._refresh_layer(layer) 374 QgsMessageLog.logMessage( 375 f'removed point {feature_id} from "{layer.name()}"', 376 MENU, level=Qgis.Info, 377 ) 378 self.iface.messageBar().pushMessage( 379 MENU, 380 'Removed last interpretation point.', 381 level=Qgis.Success, 382 duration=2, 383 ) 384 break 385 self._refresh_undo_state() 386 387 def _find_profile_canvas(self): 388 main_window = self.iface.mainWindow() 389 canvases = main_window.findChildren(QgsElevationProfileCanvas) 390 for canvas in canvases: 391 if canvas.isVisible(): 392 return canvas 393 return canvases[0] if canvases else None 394 395 def _target_layer(self, crs): 396 active = self.iface.activeLayer() 397 if _is_suitable_target(active): 398 return active 399 return self._interpretation_layer(crs) 400 401 def _search_layer(self): 402 """The layer points are added to, without creating the fallback layer.""" 403 active = self.iface.activeLayer() 404 if _is_suitable_target(active): 405 return active 406 if self._layer is not None and QgsProject.instance().mapLayer(self._layer.id()): 407 return self._layer 408 return None 409 410 def _click_to_point(self, event): 411 """(distance, elevation, xy in canvas CRS) under the click, or None off the section.""" 412 canvas = self._canvas 413 profile_point = canvas.canvasPointToPlotPoint(QPointF(event.pos())) 414 distance = profile_point.distance() 415 elevation = profile_point.elevation() 416 if math.isnan(distance) or math.isnan(elevation): 417 return None 418 419 curve = canvas.profileCurve() 420 if curve is None: 421 return None 422 point_geom = QgsGeometry(curve.clone()).interpolate(distance) 423 if point_geom.isEmpty(): 424 return None 425 return distance, elevation, point_geom.asPoint() 426 427 def _find_point_at(self, event, layer): 428 """Feature id in layer drawn within HIT_TOLERANCE_PX of the click, or None.""" 429 canvas = self._canvas 430 curve = canvas.profileCurve() 431 if curve is None: 432 return None 433 curve_geom = QgsGeometry(curve.clone()) 434 tolerance = canvas.tolerance() if hasattr(canvas, 'tolerance') else 0.0 435 candidates = [] 436 for feature in layer.getFeatures(): 437 point = feature.geometry().constGet() 438 if not isinstance(point, QgsPoint): 439 continue 440 x, y = _project_xy(point, layer.crs(), canvas.crs()) 441 xy_geom = QgsGeometry(QgsPoint(x, y)) 442 if tolerance > 0 and curve_geom.distance(xy_geom) > tolerance: 443 continue 444 distance = curve_geom.lineLocatePoint(xy_geom) 445 pixel = canvas.plotPointToCanvasPoint(QgsProfilePoint(distance, point.z())) 446 candidates.append((feature.id(), pixel.x(), pixel.y())) 447 pos = event.pos() 448 return _nearest_hit(candidates, pos.x(), pos.y(), HIT_TOLERANCE_PX) 449 450 def _on_pick(self, event): 451 if self._mode == 'edit': 452 self._edit_click(event) 453 elif self._mode == 'delete': 454 self._delete_click(event) 455 else: 456 self._add_point(event) 457 458 def _add_point(self, event): 459 picked = self._click_to_point(event) 460 if picked is None: 461 return 462 distance, elevation, xy = picked 463 canvas = self._canvas 464 465 layer = self._target_layer(canvas.crs()) 466 x, y = _project_xy(xy, canvas.crs(), layer.crs()) 467 feature = QgsFeature(layer.fields()) 468 feature.setGeometry(QgsGeometry(QgsPoint(x, y, elevation))) 469 470 feature_id = _next_feature_id(layer) 471 attrs = _point_attributes( 472 {f.name() for f in layer.fields()}, 473 feature_id, distance, elevation, 474 ) 475 for k, v in attrs.items(): 476 feature[k] = v 477 478 ok, added = layer.dataProvider().addFeatures([feature]) 479 if not ok: 480 self.iface.messageBar().pushMessage( 481 MENU, 482 'Failed to add feature to layer.', 483 level=Qgis.Warning, 484 ) 485 return 486 self._history.append((layer, _added_feature_id(added, feature))) 487 self._refresh_undo_state() 488 self._refresh_layer(layer) 489 490 id_label = f' #{feature_id}' if feature_id is not None else '' 491 QgsMessageLog.logMessage( 492 f'placed point{id_label} at x={x:.3f} y={y:.3f} z={elevation:.3f} ' 493 f'(distance={distance:.3f})', 494 MENU, level=Qgis.Info, 495 ) 496 self.iface.messageBar().pushMessage( 497 MENU, 498 f'Point at x={xy.x():.1f}, y={xy.y():.1f}, ' 499 f'elevation={elevation:.2f}', 500 level=Qgis.Success, 501 duration=2, 502 ) 503 504 def _edit_click(self, event): 505 if self._selected is None: 506 self._select_for_edit(event) 507 else: 508 self._move_selected(event) 509 510 def _select_for_edit(self, event): 511 layer = self._search_layer() 512 fid = self._find_point_at(event, layer) if layer is not None else None 513 if fid is None: 514 self.iface.messageBar().pushMessage( 515 MENU, 'No point near the click.', level=Qgis.Info, duration=2) 516 return 517 if not _can_move(layer): 518 self.iface.messageBar().pushMessage( 519 MENU, 520 f'Layer "{layer.name()}" does not support moving features.', 521 level=Qgis.Warning, 522 ) 523 return 524 layer.selectByIds([fid]) 525 self._selected = (layer, fid) 526 self.iface.messageBar().pushMessage( 527 MENU, 'Point selected; click its new position.', level=Qgis.Info, duration=2) 528 529 def _move_selected(self, event): 530 layer, fid = self._selected 531 if QgsProject.instance().mapLayer(layer.id()) is None: 532 self._selected = None 533 return 534 picked = self._click_to_point(event) 535 if picked is None: 536 return 537 distance, elevation, xy = picked 538 x, y = _project_xy(xy, self._canvas.crs(), layer.crs()) 539 provider = layer.dataProvider() 540 geometry = QgsGeometry(QgsPoint(x, y, elevation)) 541 if not provider.changeGeometryValues({fid: geometry}): 542 self.iface.messageBar().pushMessage( 543 MENU, 'Failed to move the point.', level=Qgis.Warning) 544 return 545 546 fields = layer.fields() 547 attrs = _point_attributes({f.name() for f in fields}, None, distance, elevation) 548 attrs.pop('id', None) 549 if attrs and provider.capabilities() & _CHANGE_ATTRIBUTES_CAP: 550 provider.changeAttributeValues( 551 {fid: {fields.indexOf(k): v for k, v in attrs.items()}}) 552 553 self._clear_selection() 554 self._refresh_layer(layer) 555 QgsMessageLog.logMessage( 556 f'moved point {fid} in "{layer.name()}" to x={x:.3f} y={y:.3f} ' 557 f'z={elevation:.3f} (distance={distance:.3f})', 558 MENU, level=Qgis.Info, 559 ) 560 self.iface.messageBar().pushMessage( 561 MENU, 562 f'Moved point to x={xy.x():.1f}, y={xy.y():.1f}, ' 563 f'elevation={elevation:.2f}', 564 level=Qgis.Success, 565 duration=2, 566 ) 567 568 def _delete_click(self, event): 569 layer = self._search_layer() 570 fid = self._find_point_at(event, layer) if layer is not None else None 571 if fid is None: 572 self.iface.messageBar().pushMessage( 573 MENU, 'No point near the click.', level=Qgis.Info, duration=2) 574 return 575 if not _can_delete(layer): 576 self.iface.messageBar().pushMessage( 577 MENU, 578 f'Layer "{layer.name()}" does not support deleting features.', 579 level=Qgis.Warning, 580 ) 581 return 582 if not layer.dataProvider().deleteFeatures([fid]): 583 self.iface.messageBar().pushMessage( 584 MENU, 'Failed to delete the point.', level=Qgis.Warning) 585 return 586 self._history = [ 587 (lyr, f) for lyr, f in self._history 588 if not (lyr.id() == layer.id() and f == fid) 589 ] 590 self._refresh_undo_state() 591 self._refresh_layer(layer) 592 QgsMessageLog.logMessage( 593 f'deleted point {fid} from "{layer.name()}"', MENU, level=Qgis.Info) 594 self.iface.messageBar().pushMessage( 595 MENU, 'Deleted interpretation point.', level=Qgis.Success, duration=2) 596 597 def _interpretation_layer(self, crs): 598 project = QgsProject.instance() 599 if self._layer is not None and project.mapLayer(self._layer.id()): 600 return self._layer 601 602 uri = f'PointZ?crs={crs.authid()}' if crs.isValid() else 'PointZ' 603 layer = QgsVectorLayer(uri, LAYER_NAME, 'memory') 604 provider = layer.dataProvider() 605 provider.addAttributes([ 606 QgsField('id', QMetaType.Type.Int), 607 QgsField('distance', QMetaType.Type.Double), 608 QgsField('elevation', QMetaType.Type.Double), 609 ]) 610 layer.updateFields() 611 project.addMapLayer(layer) 612 self._layer = layer 613 return layer