Obtenez toutes les lignes qui entourent un point

J'utilise QGIS et j'ai un point et le réseau routier. J'ai besoin d'extraire automatiquement les noms de routes qui entourent le point spécifique. L'analyse du voisin le plus proche et les zones tampons ne peuvent pas faire le travail, car dans de nombreux cas, le point est plus proche, en termes de distance mesurée, des routes voisines et non des routes environnantes. Y a-t-il des idées sur la façon dont seules les routes environnantes pourraient être extraites?

À propos de mes données de test:

1. Comme OSM Road Data, chaque géométrie de route se termine à un carrefour.
2. Chaque route a un identifiant unique.

SOLUTION I

S'il y a les deux hypothèses:

1. Les routes construisent des quartiers.

2. Vous travaillez dans un système métrique.

L'idée est d'augmenter / diminuer les coordonnées X et Y du point. Si vous travaillez dans un système métrique, vous pouvez vous rendre à 1 m à l'est de votre point, créer un nouveau point et créer une ligne avec le point d'origine. Vous allez plus loin vers l'est jusqu'à ce que la ligne croise une route. Pour rechercher une intersection à l'ouest, vous devez soustraire 1 m de la coordonnée X d'origine. De même pour la coordonnée Y. S'il n'y a pas de route au Nord / Est / Sud / Ouest le compteur s'arrête à 1000 (m). Lorsque vous savez qu'il pourrait y avoir une route à une distance supérieure à 1000 mètres, vous devez modifier cette valeur.

Vous pouvez résoudre la tâche avec le code suivant:

Édité

for lyr in QgsMapLayerRegistry.instance().mapLayers().values():
    if lyr.name() == "point":
        startpoint = lyr

for lyr in QgsMapLayerRegistry.instance().mapLayers().values():
    if lyr.name() == "roads":
        roads = lyr

startpoint_iter = startpoint.getFeatures()
for feature in startpoint_iter:
    geom = feature.geometry()
    if geom.type() == QGis.Point:
        xy = geom.asPoint()
        x,y = xy[0], xy[1]

line_start = QgsPoint(x,y)      

def reached(direction, count_m):
    road_reached = None
    road = None
    count=1
    while road_reached < 1 and count <=count_m:
        count += 1
        if direction == 'N':
            line_end = QgsPoint(x, y+count)
        if direction == 'E':
            line_end = QgsPoint(x+count,y)
        if direction == 'S':
            line_end = QgsPoint(x,y-count)
        if direction == 'W':
            line_end = QgsPoint(x-count,y)
        line = QgsGeometry.fromPolyline([line_start,line_end])
        for f in roads.getFeatures():
            if line.intersects(f.geometry()):
                road_reached = 1
                road = f['name']
                print road

reached('N', 1000)
reached('E', 1000)
reached('S', 1000)
reached('W', 1000)

Un autre exemple pour montrer que la route e à l'Est n'est pas reconnue comme une route voisine du point.

Comment appeler la fonction et la sortie:

>>>>reached('N', 1000)
road a
>>>>reached('E', 1000)
road b
>>>>reached('S', 1000)
road c
>>>>reached('W', 1000)
road d

S'il y a plus de 4 routes entourant le point, vous devez regarder dans plus de directions (changer à la fois X et Y). Vous pouvez également modifier l'azimut de votre ligne, ce qui signifie que vous pouvez la faire pivoter d'un degré dans la plage de 0 à 360 °.

SOLUTION II

Inspiré du commentaire, vous pouvez également Polygonizecommencer par vos routes. Pour cela vous pouvez utiliser un outil de QGIS: Processing > Toolbox > QGIS geoalgorithms > Vector geometry tools > Polygonize. Renommez le calque temporaire en polygon. En supposant que vous souhaitez uniquement avoir les noms de route pour le point qui est complètement entouré de routes. Sinon , vous devez utiliser SOLUTION I . Cela ne fonctionne que si toutes les routes sont connectées (cassées)!

Le point doit d'abord intersecter avec le polygone. L'idée est maintenant que les deux ANDpoints de départ d'une ligne englobante doivent se croiser avec le polygone.

for lyr in QgsMapLayerRegistry.instance().mapLayers().values():
    if lyr.name() == "point":
        startpoint = lyr

for lyr in QgsMapLayerRegistry.instance().mapLayers().values():
    if lyr.name() == "polygon":
        poly = lyr

for lyr in QgsMapLayerRegistry.instance().mapLayers().values():
    if lyr.name() == "roads":
        roads = lyr

for h in startpoint.getFeatures():
    for g in poly.getFeatures():
        if h.geometry().intersects(g.geometry()):
            poly_geom = g.geometry()
            for f in roads.getFeatures():
                geom = f.geometry().asPolyline()
                start_point = QgsGeometry.fromPoint(QgsPoint(geom[0]))
                end_point = QgsGeometry.fromPoint(QgsPoint(geom[-1]))
                if poly_geom.intersects(start_point) and poly_geom.intersects(end_point):
                    print f['name']

Le résultat:

road c
road b
road e
road f