Rotating point objects toward lines using QGIS

Question

I have a set of points that I'd like to have rotated toward nearest line by adding an attribute column angle.

Original points with no rotation

Points rotated toward line object

I have found and tweaked this script to somewhat achieve what I want. The only problem is that this method only does a search using a rectangular box which sometimes gives me a less satisfying result. Is there a way to do the search using a circle instead?

Here's my tweaked script:

from math import atan2
get layers
buildings = QgsProject.instance().mapLayersByName('buildings')[0]
roads = QgsProject.instance().mapLayersByName('roads')[0]
create angle field
fieldname = "angle"
buildings.dataProvider().addAttributes([QgsField(fieldname,QVariant.Int)])
buildings.updateFields()
get fieldindex
fni = buildings.fields().indexFromName('angle')
initializations
tolerance = 30  # search tolerance
buildings.startEditing()
loop over all buildings
for building in buildings.getFeatures():
    x = building.geometry().asPoint().x()
    y = building.geometry().asPoint().y()
# get the rectangular search area 
searchRect = QgsRectangle(x - tolerance, y - tolerance,  x + tolerance, y + tolerance)

# find roads 
for road in roads.getFeatures(QgsFeatureRequest().setFilterRect(searchRect)):
    # get the nearest vertex on road and the one before and after
    pnt, v, b, a, d = road.geometry().closestVertex(building.geometry().asPoint())
    p1 = road.geometry().vertexAt(v)
    # when vertex before exists look back, otherwise look forward
    if v>-1 and b>-1:
        p2 = road.geometry().vertexAt(b)
    elif v>-1 and a>-1:
        p2 = road.geometry().vertexAt(a)

    # calculate azimuth
    angle = atan2(p2.x() - p1.x(), p2.y() - p1.y()) / 0.017453
    building[fni] = angle
    buildings.updateFeature(building)


save changes and stop editing
buildings.commitChanges()

Answer 1

I'm not sure if this will give you more satisfactory results than what you are currently using, but I had a play around with some test data and came up with the following script which you are welcome to try. Instead of using a search rectangle or radius, it simply calculates the rotation from the azimuth of the nearest segment of the nearest line feature for each point feature.

# get layers
buildings = QgsProject.instance().mapLayersByName('buildings')[0]
roads = QgsProject.instance().mapLayersByName('roads')[0]

# create angle field
fieldname = "angle"
buildings.dataProvider().addAttributes([QgsField(fieldname,QVariant.Int)])
buildings.updateFields()

# get fieldindex
fni = buildings.fields().indexFromName('angle')

# build spatial index on roads layer
index = QgsSpatialIndex(roads.getFeatures(), flags=QgsSpatialIndex.FlagStoreFeatureGeometries)

att_map = {}

# loop over all buildings
for building in buildings.getFeatures():
    # find nearest road feature to current building
    closest_feature_ids = index.nearestNeighbor(building.geometry())
    if closest_feature_ids:
        if len(closest_feature_ids) == 1:
            closest_feature = roads.getFeature(closest_feature_ids[0])
        else:
            nearest_points = [f.geometry().nearestPoint(building.geometry()) for f in roads.getFeatures(closest_feature_ids)]
            closest = [QgsGeometry.collectGeometry(nearest_points).nearestPoint(building.geometry())][0]
            closest_feature = [f for f in roads.getFeatures(closest_feature_ids) if f.geometry().intersects(closest.buffer(0.001, 8))][0]
    # find closest line segment of road feature
    a, b, c, d = closest_feature.geometry().closestSegmentWithContext(building.geometry().asPoint())
    # get azimuth of line from vertices before and after closest segment
    azimuth = closest_feature.geometry().vertexAt(c).azimuth(closest_feature.geometry().vertexAt(c-1))
    # create attribute map of feature ids, field index and azimuth value
    att_map[building.id()] = {fni: azimuth}

# update attribute values in buildings layer
buildings.dataProvider().changeAttributeValues(att_map)

Before:

After (with marker rotation using data-defined values from 'angle' field):

You can also try this slightly different script which does not use a spatial index so it will be slower on large datasets but may give a more satisfactory result (I don't know).

# get layers
buildings = QgsProject.instance().mapLayersByName('buildings')[0]
roads = QgsProject.instance().mapLayersByName('roads')[0]
create angle field
fieldname = "angle"
buildings.dataProvider().addAttributes([QgsField(fieldname,QVariant.Int)])
buildings.updateFields()
get fieldindex
fni = buildings.fields().indexFromName('angle')
att_map = {}
loop over all buildings
for building in buildings.getFeatures():
    # find nearest road feature to current building
    nearest_points = [f.geometry().nearestPoint(building.geometry()) for f in roads.getFeatures()]
    closest = [QgsGeometry.collectGeometry(nearest_points).nearestPoint(building.geometry())][0]
    closest_feature = [f for f in roads.getFeatures() if f.geometry().intersects(closest.buffer(0.001, 8))][0]
    # find closest line segment of road feature
    a, b, c, d = closest_feature.geometry().closestSegmentWithContext(building.geometry().asPoint())
    # get azimuth of line from vertices before and after closest segment
    azimuth = closest_feature.geometry().vertexAt(c).azimuth(closest_feature.geometry().vertexAt(c-1))
    # create attribute map of feature ids, field index and azimuth value
    att_map[building.id()] = {fni: azimuth}
update attribute values in buildings layer
buildings.dataProvider().changeAttributeValues(att_map)