2

I wrote a such script to smooth my line. It checks whether a line of two points intersects the border of a polygon and, if not, then draws a line.

This script works well for a simple example, but bad for this example. For some reason it always reduces the line to two points (start and end) and the intersection condition is not met. The boundary line is a Multiline in both cases.

Can anyone help me fix the condition so that the script runs correctly?

point_layer = QgsProject.instance().mapLayersByName("888 Vertices")[0]
path=[]
for feature in point_layer.getFeatures():
    geom = feature.geometry()
    point = geom.asPoint()
    path.append((point.x(), point.y()))
    print(point)
print(len(path))


bound_layer = QgsProject.instance().mapLayersByName("9 boundary 3 collected")[0]
feature = bound_layer.getFeature(1)

if len(path) > 2:
    k=0
    t=[]
    t.append(path[0])
    for i in range(1, len(path)-1):   
        if QgsGeometry.fromPolylineXY([QgsPointXY(t[k][0], t[k][1]),QgsPointXY(path[i+1][0], path[i+1][1])]).intersects(feature.geometry()):
            k+=1
            t.append(path[i])
    k+=1
    t.append(path[-1])
    
print(len(t))    
print(t)


vl = QgsVectorLayer(f"LineString?crs={point_layer.crs().authid()}", "smoothed_path", "memory")

pr = vl.dataProvider()

# Create a new feature
seg = QgsFeature()

# Add the shortest path to the feature as a polyline
seg.setGeometry(QgsGeometry.fromPolyline([QgsPoint(pt[0], pt[1]) for pt in t]))

# Add the feature to the layer
pr.addFeatures([seg])

# Update extent of the layer
vl.updateExtents()
QgsProject.instance().addMapLayer(vl)

For this example it works good and really reduces the source line:

enter image description here

But it doesn't work for this example at all (orange = boundary line, yellow = source line and blue = smoothed line, it's incorrect). Blue line should be at least the same as yellow line but a little but smoothed, without unnecessary vertices.

enter image description here

I have loaded my data in here: https://drive.google.com/drive/folders/1NH4wE1u8wii5c1MK-s17deC-QLDj1Ox9?usp=sharing

0

1 Answer 1

1

There are several suggestions regarding your code:

  • Break your '9 boundary 3 collected' lines with the "Multipart to singleparts":

    bound_layer = QgsProject.instance().mapLayersByName("9 boundary 3 collected")[0]
    params = {
       'INPUT': bound_layer,
       'OUTPUT': 'TEMPORARY_OUTPUT'
       }
    bound_layer_single = processing.run("native:multiparttosingleparts", params)['OUTPUT']
    bound_features = [feature.geometry() for feature in bound_layer_single.getFeatures()]
    
  • check each feature for intersection separately by means of the Python's built-in function any()

    if any(line.intersects(bound_feature) for bound_feature in bound_features):
    
  • append values to the path = [] directly as a QgsPointXY-object, no need to create a tuple with coordinates i.e. (point.x(), point.y()) and after reconstructing them back. Then your code will be also simplified in this line:

    if QgsGeometry.fromPolylineXY([QgsPointXY(t[k][0], t[k][1]), path[i+1]]).intersects(feature.geometry())
    
  • Python - List Comprehension is a great asset to create the path-object:

    path = [QgsPointXY(feature.geometry().asPoint()) for feature in point_layer.getFeatures()]
    print(len(path)) #58
    
  • The t-list with its QgsPoint-objects can be adjusted via the map() function combined with the lambda:

    points = list(map(lambda point: QgsPoint(point), t))
    # Add the shortest path to the feature as a polyline
    seg.setGeometry(QgsGeometry.fromPolyline(points))
    
  • Sometimes the order of coordinates of a point can swap, therefore use the corresponding methods e.g. x() and y() of the QgsPointXY class:

    if QgsGeometry.fromPolylineXY([QgsPointXY(t[k].x(), t[k].y()), path[i+1]]).intersects(feature.geometry()):
    

    But the QgsGeometry-object you intended to achieve can be initiated even easier:

    if QgsGeometry.fromPolylineXY([t[k], path[i+1]]).intersects(feature.geometry()):
    

    See, that the t[k]-variable is of the QgsPointXY-class type. The QgsGeometry::fromPolylineXY creates a new LineString geometry from a list of QgsPointXY points.

  • Imports can be useful to make your code reusable in the future in standalone applications:

    # imports
    import processing
    from qgis.core import (QgsProject, QgsVectorLayer, QgsPointXY, QgsPoint, QgsFeature, QgsGeometry)
    
  • leave comments in your code

So, in the end, your code may look like this:

# imports
import processing
from qgis.core import (QgsProject, QgsVectorLayer, QgsPointXY, QgsPoint, QgsFeature, QgsGeometry)

point_layer = QgsProject.instance().mapLayersByName("888 Vertices")[0]

path = [QgsPointXY(feature.geometry().asPoint()) for feature in point_layer.getFeatures()]

bound_layer = QgsProject.instance().mapLayersByName("9 boundary 3 collected")[0]
params = {
    'INPUT': bound_layer,
    'OUTPUT': 'TEMPORARY_OUTPUT'
    }
bound_layer_single = processing.run("native:multiparttosingleparts", params)['OUTPUT']
bound_features = [feature.geometry() for feature in bound_layer_single.getFeatures()]

if len(path) > 2:
    k, t = 0, []
    t.append(path[0])
    for i in range(1, len(path) - 1):
        line = QgsGeometry.fromPolylineXY([t[k], path[i + 1]])
        if any(line.intersects(bound_feature) for bound_feature in bound_features):
            k += 1
            t.append(path[i])
    k += 1
    t.append(path[-1])

vl = QgsVectorLayer(f"LineString?crs={point_layer.crs().authid()}", "smoothed_path", "memory")

pr = vl.dataProvider()

# Create a new feature
seg = QgsFeature()

points = list(map(lambda point: QgsPoint(point), t))
# Add the shortest path to the feature as a polyline
seg.setGeometry(QgsGeometry.fromPolyline(points))

# Add the feature to the layer
pr.addFeatures([seg])

# Update extent of the layer
vl.updateExtents()
QgsProject.instance().addMapLayer(vl)

And then it should work as supposed.

result

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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