1

I have two layers: point and line (postgis or shp) The line layer represent a hiking trail. The point layer the mountain passes. The points are exactly on the same position like the node of the line. Now I would like to rotate the points in orientation of the line. Is there a way to do that automatically?

enter image description here

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
MartinMap
  • 8,262
  • 7
  • 50
  • 114

1 Answers1

8

You need to determine the azimuth of the polyline segment next to the point. It is possible to calculate this value and save it to an attribute, or calculate it as an expression which give the azimuth as result.

How to use the following Python script:

  1. Point layer and line layer are added to the legend (in the script named bridges and trails) and change the names in the code
  2. Add an attribut of type REAL to the point layer (in the script named angle), when using a different name, change the code in line 8
  3. Open Python Console -> Editor and copy & paste the code
  4. Run the script
from math import atan2

# get layers
bridges = QgsMapLayerRegistry.instance().mapLayersByName('bridges')[0]
trails = QgsMapLayerRegistry.instance().mapLayersByName('trails')[0]

# get fieldindex
fni = bridges.fieldNameIndex('angle')

# initializations
tolerance = 0.01  # search tolerance
bridges.startEditing()

# loop over all bridges
for bridge in bridges.getFeatures():
    x = bridge.geometry().asPoint().x()
    y = bridge.geometry().asPoint().y()

    # get the rectangular search area 
    searchRect = QgsRectangle(x - tolerance, y - tolerance,  x + tolerance, y + tolerance)

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

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

# save changes and stop editing
bridges.commitChanges()

The value for variable tolerance depends on the units of length measurements. Tt may be necessary to choose the appropriate value to really pick the nearest vertex.

Use this attribut angle as variable for data defined rotation of the symbols:

rotated symbols

You may wish to complicate things a little bit and calculate azimuths to the vertex before AND after and take the average.

The dynamic expression version:

from qgis.core import *
from qgis.gui import *
from math import atan2

@qgsfunction(args='auto', group='Custom')
def get_azimuth(layername, tolerance, feature, parent):
    # initializations
    trails  = QgsMapLayerRegistry.instance().mapLayersByName(layername)[0]
    x = feature.geometry().asPoint().x()
    y = feature.geometry().asPoint().y()

    # get the rectangular search area 
    searchRect = QgsRectangle(x - tolerance, y - tolerance,  x + tolerance, y + tolerance)

    # find trails 
    for trail in trails.getFeatures(QgsFeatureRequest().setFilterRect(searchRect)):
        # get the nearest vertex on trail and the one before and after
        pnt, v, b, a, d = trail.geometry().closestVertex(feature.geometry().asPoint())
        if v == -1:
            continue
        p1 = trail.geometry().vertexAt(v)
        # when vertex before exists look back, otherwise look forward
        if v>-1 and b>-1:
            p2 = trail.geometry().vertexAt(b)
        elif v>-1 and a>-1:
            p2 = trail.geometry().vertexAt(a)
        # calculate azimuth
        angle = atan2(p2.x() - p1.x(), p2.y() - p1.y()) / 0.017453
        return angle

expression

field calculator

Detlev
  • 4,608
  • 19
  • 26
  • Looks pretty nice! But I'm sorry, I don't know how to run that script in qgis. I would like to get the result directly as expression in the Rotation-Filed. – MartinMap Jan 29 '16 at 19:18
  • 1
    @MAP I added the answer, hope you get it. Sometimes I had to wait several seconds for the symbols to get correct rotated. It named the layer bridges but you seem to be interested in passes, don't you? But the principle remains the same. – Detlev Jan 29 '16 at 20:21
  • Great, works like it should! You've saved me a lot of work! Thanks! – MartinMap Jan 30 '16 at 11:47
  • I've used this technique to build an expression, but it doesn't quite work correctly. The cross lines (point symbols) should be perpendicular to the stream, and some of them appear correct while others are closer to parallel. Any idea where my problem is?

    https://i.stack.imgur.com/0r8Wv.png

    – Nathan Oct 30 '17 at 04:15
  • @Nathan If you could share the data example from the picture then I would use this to examine the problem. As I wrote in the answer I encountered this problem in very few cases, but you have many of them. So chances are good to find a solution. – Detlev Oct 30 '17 at 11:02
  • The data source is the U.S. National Hydrography Dataset, and I'm using an extract for New York State: http://prd-tnm.s3-website-us-west-2.amazonaws.com/?prefix=StagedProducts/Hydrography/NHD/State/HighResolution/GDB/ Filename: NHD_H_New_York_GDB.zip

    In my example, I've opened the GDB and loaded the Flowline layer for my streams, and the Points layer to get rapids and falls along the stream. (I've filtered the Points layer to "FType" IN (431,487) to get just those features.)

    – Nathan Oct 30 '17 at 12:37
  • To find this particular cluster of points, use extent coords 581507,4984943 : 585626,4987593. – Nathan Oct 30 '17 at 12:41
  • 1
    The variable tolerance is used to search lines near to the points. Choose this value as small as possible, otherwise the script has to be changed to evaluate for the closest vertex in the search rectangle. With the data you mentioned (in NAD83) I get reasonable results with tolerance equal to 0.0001. Keep in mind that this value depends on the projection respectively length units used. – Detlev Oct 30 '17 at 17:17
  • Aha—tolerance was indeed the key! So, if I understand correctly: I had thought 0.01 would be more than sufficient, since my project is in meters (UTM) and so 0.01 equals a centimeter. But the hydrography layer is in unprojected NAD 83, so its units are degrees—and 0.0001 equals about half a second. So I needed to think about the units of the layer, not of the project overall. – Nathan Nov 01 '17 at 15:14