Cut a linear feature according to the directions or azimuths

Question

I have a shapefile linear that contains a single items as in the image and I want to divide it into several section respecting the azimuth and the line direction.

i have this code :

    import os
    import fiona
    from shapely.geometry import mapping
    from shapely.geometry import Point,LineString
    import math

    os.chdir('SHP')

    def azimut(point1, point2):
           '''Retourne l'azimuth de la ligne entre 2 points shapely'''
           angle = math.atan2(point2.x - point1.x, point2.y - point1.y)
           return math.degrees(angle)
           print angle
           **if angle &amp;&gt;= 0: else math.degrees(angle) + 180** **Error syntax**


    **def paires(liste):
         '''parcourt une liste par paires'''
         for i in range(1, len(liste)):
             yield liste[i-1], liste[i]**

    with fiona.collection('BOM1.shp', 'r') as entree:
        # copie du schema de la couche et création d'un nouveau champ 'azimut'
        schema = entree.schema.copy()
        schema['properties']['azimut'] = 'int'
        # création d' une nouvelle couche avec le schéma résultant
        with fiona.collection('testligne_azim.shp', 'w', 'ESRI Shapefile', schema) as sortie:
            for ligne in entree:
            # utilisation de la fonction paire() pour extraire les segments de lignes
                for seg_start, seg_end in paires(ligne['geometry']['coordinates']):
                    # création d'une ligne en fonction des points des segments
                    line_start =Point(seg_start)
                    line_end = Point(seg_end)
                    segment = LineString([line_start.coords[0],line_end.coords[0]])
                    # copie des attributs d'entrée et ajout de la valeur résultante
                    # de la fonction azimut()
                    elem = {}                  
                    elem['properties'] = ligne['properties'] 
                    elem['properties']['azimut'] = azimut(line_start, line_end)
                    elem['geometry'] = mapping(segment)
                    sortie.write(elem)

but the problem is in :

if angle &amp;&gt;= 0: else math.degrees(angle) + 180

i have syntax error angle &>= 0

and in:

def paires(liste):
         '''parcourt une liste par paires'''
         for i in range(1, len(liste)):
             yield liste[i-1], liste[i]

i have this error position 29: ordinal not in range(128)

why &>= instead of >= and does angle contains non-ascii char (e.g. degree symbol)? — Learner, Feb 10 '16 at 05:09

score 1 · Answer 1 · edited Apr 13 '17 at 12:33

It is a simple transliteration (HTML) error ("&>" = > in HTML).
The script comes from Python: traitement des couches vectorielles dans une perspective géologique, lecture et enregistrement des couches sous forme de dictionnaires avec le module Fiona and the complete code in English is given in How do I find vector line bearing in QGIS or GRASS?

 def azimuth(point1, point2):
   '''azimuth between 2 shapely points (interval 0 - 180°, for my work)'''
   angle = math.atan2(point2.x - point1.x, point2.y - point1.y)
   return math.degrees(angle)if angle>0 else math.degrees(angle) + 180

def pair(list):
   '''Iterate over pairs in a list '''
    for i in range(1, len(list)):
        yield list[i-1], list[i]

The second function (def pair(list)) is a simple Python generator that is not the cause of your error.

Note

Honestly, knowing Python and a few angular rules, it was not too difficult to correct the first function...

You can use return math.degrees(angle) % 360

Cut a linear feature according to the directions or azimuths

1 Answers1