2

I have a vector layer in QGIS. I'm looking to return the xyz, not of the layer's feature(s), but of the intersection point between a "Select Features by Freehand" geometry and the feature(s), then return it into the attribute table. The idea would be that when the select tool's line crosses any line in the vector layer, the features gets selected and the interpolated xyz at those intersection points gets returned as three numeric fields.

Is this sort of thing possible using Field Calculator, or Actions?

I have had a play around and I can't find a way to capture the selection geometry to use in something like ST_Intersection. I've also considered a custom search tool, but it seems a bit tricky.

I haven't created a custom plugin or coded a function yet, but this would definitely be an option if there's a solution there, but how would I get the selection tool's geometry or is this only possible by creating a new layer which holds the query geometry?

2

Interesting question. In my opinion, a Python script which implements a custom QgsMapTool class is the way to go here. I don't know another way to interact with the geometry of a QgsRubberBand drawn on the map canvas.

I have developed a script which (I think) should do what you want which you are welcome to use and adapt. It could also be turned into a small plugin.

Please note the following points:

  • Your question is not entirely explicit about which geometry type your vector layer has, but I am assuming lines. If you have polygons with z values stored in the vertices, the script will need some modification.

  • The map tool works the same as the native freehand selection tool; click once to start digitizing, click again to finish.

  • Since you mentioned one new field for x, y & z intercept values, the script is designed to work when the digitized polygon crosses a line feature exactly once. If you do cross a feature more than once, only one of the intersection points will be saved (I think it's the one closest to the start point of the linestring). To help you know which intersection point's XYZ values are saved to the attribute table, each saved point is visualized with a QgsVertexMarker when you finish digitizing the rubber band.

I tested this on a line layer with MultiLineStringZ geometry:

enter image description here

Please see below for a gif showing how the tool works:

enter image description here

The full script is below. To use it, open the Python console using the enter image description here button (or Ctrl+Alt+P) and open a new editor with this button:

enter image description here

Paste in the script and click Run:

enter image description here

  • Make sure your line layer is selected as the active layer and start digitizing on the canvas with the freehand tool.
#-----------------------------------------------------------
# Copyright (C) 2021 Ben Wirf
#-----------------------------------------------------------
# Licensed under the terms of GNU GPL 2
#
# 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.
#---------------------------------------------------------------------

class DrawFreehandPolygon(QgsMapTool):
    
    def __init__(self, canvas):
        self.canvas = canvas
        QgsMapTool.__init__(self, self.canvas)
        self.polygon_rb = QgsRubberBand(self.canvas, QgsWkbTypes.PolygonGeometry)
        self.set_points = []
        self.vertex_markers = []
        self.drawing_mode = False
        
        
    def canvasPressEvent(self, e):
        self.set_points.append(e.mapPoint())
        if self.drawing_mode == False:
            self.clear_markers()
            self.polygon_rb.reset()
            self.drawing_mode = True
        elif self.drawing_mode:
            self.clear_markers()
            self.polygon_rb.reset()
            pts_arr = [p for p in self.set_points]
            self.polygon_rb.setToGeometry(QgsGeometry.fromPolygonXY([pts_arr]))
            self.polygon_rb.setStrokeColor(QColor(240, 165, 105))
            self.polygon_rb.setFillColor(QColor(240, 165, 105, 85))
            self.polygon_rb.setWidth(2)
            self.polygon_rb.show()
            
            x_geom = self.transform_polygon(QgsProject.instance().crs(),
                                iface.activeLayer().crs(),
                                QgsGeometry.fromPolygonXY([pts_arr]))
            
            self.intersection_points(x_geom)
            self.drawing_mode = False
            self.set_points.clear()

        
    def canvasMoveEvent(self, e):
        if self.drawing_mode:
            self.polygon_rb.reset()
            self.set_points.append(e.mapPoint())
            if len(self.set_points) > 2:
                pts_arr = [p for p in self.set_points]
                self.polygon_rb.setToGeometry(QgsGeometry.fromPolygonXY([pts_arr]))
                self.polygon_rb.setStrokeColor(QColor(240, 165, 105))
                self.polygon_rb.setFillColor(QColor(240, 165, 105, 85))
                self.polygon_rb.setWidth(2)
                self.polygon_rb.show()

    def transform_polygon(self, canvas_crs, layer_crs, g):
        if canvas_crs != layer_crs:
            xform = QgsCoordinateTransform(canvas_crs, layer_crs, QgsProject.instance())
            g.transform(xform)
        return g

    def intersection_points(self, rb_geom):
        lyr = iface.activeLayer()
        current_fld_names = [fld.name() for fld in lyr.fields()]
        flds_to_add = []
        if not 'X' in current_fld_names:
            flds_to_add.append(QgsField('X', QVariant.Double, '', len=25, prec=6))
        if not 'Y' in current_fld_names:
            flds_to_add.append(QgsField('Y', QVariant.Double, '', len=25, prec=6))
        if not 'Z' in current_fld_names:
            flds_to_add.append(QgsField('Z', QVariant.Double, '', len=25, prec=3))
        if flds_to_add:
            lyr.dataProvider().addAttributes(flds_to_add)
            lyr.updateFields()
        x_idx = lyr.fields().lookupField('X')
        y_idx = lyr.fields().lookupField('Y')
        z_idx = lyr.fields().lookupField('Z')
        ### Clear existing X,Y,Z attributes
        for f in lyr.getFeatures():
            lyr.dataProvider().changeAttributeValues({f.id(): {x_idx: NULL, y_idx: NULL, z_idx: NULL}})
        ###
        intersections = []
        for f in lyr.getFeatures():
            if f.geometry().intersects(rb_geom):
                intersections.append(f.geometry().intersection(rb_geom))
        for geom in intersections:
            feat_to_update = [f for f in lyr.getFeatures() if f.geometry().intersects(geom)][0]
            vert_pnts = [v for v in geom.vertices() if not QgsGeometry().fromPointXY(QgsPointXY(v)).buffer(0.0000001, 8).within(rb_geom)]
            if vert_pnts:
                self.create_vertex_marker(vert_pnts[0])
                x_val = vert_pnts[0].x()
                y_val = vert_pnts[0].y()
                z_val = vert_pnts[0].z()
                lyr.dataProvider().changeAttributeValues({feat_to_update.id(): {x_idx: x_val, y_idx: y_val, z_idx: z_val}})
    
    def create_vertex_marker(self, pt):
        m = QgsVertexMarker(self.canvas)
        m.setCenter(QgsPointXY(pt))
        m.setColor(QColor(0, 0, 255))
        m.setIconSize(7)
        m.setIconType(QgsVertexMarker.ICON_BOX) # or ICON_CROSS, ICON_CIRCLE
        m.setPenWidth(3)
        self.vertex_markers.append(m)
    
    def clear_markers(self):
        if self.vertex_markers:
            for m in self.vertex_markers:
                self.canvas.scene().removeItem(m)
            self.vertex_markers.clear()
    
    def deactivate(self):
        self.polygon_rb.reset()
        self.clear_markers()
        

mt = DrawFreehandPolygon(iface.mapCanvas())
iface.mapCanvas().setMapTool(mt)
6
  • This is great, thank you so much for this really detailed answer!! I've pretty much managed to get this working, albeit with changing the rb type to line (and learning a lot about QGIS in the process!). I encountered a problem with getting the intersection point to re-intersect with the feature - was this why you used the buffer in your case? Also, don't worry about this too much, as can always post a new question, but I've noticed that I actually have to do the selection twice (I'm viewing "show selected features") for the attributes to update. Any idea why this might be?
    – jraclarke
    Jul 15 '21 at 16:35
  • Hi @jraclarke, can you explain why you changed the rubber band type to line? I wrote this script specifically to answer your question and I used a polygon rb because I thought you wanted something like the Select by Freehand tool. I really expected that you could paste this code into the console and it would work for you. Did you try that and how did it not suit you? Yes, you certainly could use a line rubber band, but that would completely change the logic for the rest of the script to get the intersection points.
    – Ben W
    Jul 16 '21 at 6:24
  • Wait, are you actually activating the native Select by Freehand tool? Because I'm not sure what you mean by "viewing show selected features". This script uses a custom subclass which inherits from QgsMapTool and overrides the CanvasPressEvent method to implement the custom functionality of getting the intersection points and updating the attribute table instead of selecting. It would help me help you if could provide a sample of the layer you are working on via a Google Drive or Dropbox link or similar if at all possible.
    – Ben W
    Jul 16 '21 at 6:35
  • Hi @Ben W, only changed to line because someone asked if we could also have it so that we provided a cross-section across the features as well (sorry should have specified); your code is still valid as I'm using both custom tools (cross-section and freehand)! I'm not using the native Freehand tool, but wasn't sure why they weren't "selected" in the custom, but what you mentioned in your 2nd comment makes sense of this behaviour. To be honest, this isn't a problem but was just interested. Thanks again for the massive amount of help on this, it's really appreciated!
    – jraclarke
    Jul 16 '21 at 9:45
  • 1
    Thanks @Ben W, yeah I'll add that functionality in. Line custom tool is working as well, so all good. No worries, answer accepted, and will upvote as well when I have enough reputation!
    – jraclarke
    Jul 16 '21 at 10:32
1

You can create a new polygon layer instead of the selection tool and than create the intersection points with the lines using QGIS expressions on the line layer:

intersection (
    $geometry,
    boundary (
        geometry (
            get_feature_by_id (
                'polygon',
                1
))))

In the field calculator, use x(), y() and z() and insert the expression from above to get the respective values.

Screenshot showing the expression in action, creating the intersection points with geometry generator: enter image description here

If the same line crosses the selection polygon several times, you will get only one value. Either use explode lines or Split with lines to get separate lines for each intersection.

Another possibility is to use the expression above with Geometry by expression to create the intersection points as actual geometries, than run Multipart to singleparts. On this new point-layer, use $x and $y for the coordinate values.

1
  • Great stuff and thanks. For the time being I'm doing this via selection - but some future functionality is going to involve doing it by layer, so I'll use this.
    – jraclarke
    Jul 15 '21 at 16:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.