Uzyskaj wszystkie linie, które zawierają punkt

Używam QGIS i mam sieć punktów i dróg. Muszę automatycznie wyodrębnić nazwy dróg, które zawierają konkretny punkt. Analiza najbliższych sąsiadów i stref buforowych nie może wykonać zadania, ponieważ w wielu przypadkach punkt jest bliżej, ze względu na zmierzoną odległość, do dróg sąsiadujących, a nie otaczających. Czy są jakieś pomysły, w jaki sposób można wydobyć tylko otaczające drogi?

O moich danych testowych:

1. Podobnie jak dane drogowe OSM, każda geometria drogi kończy się na rozdrożu.
2. Każda droga ma unikalny identyfikator.

ROZWIĄZANIE I

Jeśli są dwa założenia:

1. Drogi budują dzielnice.

2. Pracujesz w systemie metrycznym.

Chodzi o zwiększenie / zmniejszenie współrzędnych X i Y punktu. Jeśli pracujesz w systemie metrycznym, możesz przejść 1 m na wschód od swojego punktu, utworzyć nowy punkt i utworzyć linię z punktem oryginalnym. Jedziesz dalej na wschód, aż linia przecina drogę. Aby wyszukać skrzyżowanie na zachodzie, odejmij 1m od oryginalnej współrzędnej X. To samo dla współrzędnej Y. Jeśli na północy / wschodzie / południu / zachodzie nie ma drogi, licznik zatrzymuje się na 1000 (m). Jeśli wiesz, że droga może znajdować się w odległości większej niż 1000 m, musisz zmienić tę wartość.

Możesz rozwiązać zadanie za pomocą następującego kodu:

Edytowane

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)

Kolejny przykład pokazujący, że droga e na Wschodzie nie jest uznawana za pobliską drogę punktu.

Jak wywołać funkcję i wyjście:

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

Jeśli punkt obejmuje więcej niż 4 drogi, musisz spojrzeć w większą liczbę kierunków (zmień X i Y). Lub możesz zmienić azymut swojej linii, co oznacza, że ​​możesz obrócić ją o jeden stopień w zakresie 0-360 °.

ROZWIĄZANIE II

Zainspirowany komentarzem, możesz Polygonizenajpierw swoje drogi. Do nich można użyć narzędzia z QGIS: Processing > Toolbox > QGIS geoalgorithms > Vector geometry tools > Polygonize. Zmień nazwę warstwy tymczasowej na polygon. Zakładając, że chcesz mieć tylko nazwy dróg dla punktu, który jest całkowicie otoczony drogami. W przeciwnym razie trzeba użyć roztworu I . Działa to tylko wtedy, gdy wszystkie drogi są połączone (zatrzaśnięte)!

Najpierw punkt musi przecinać się z wielokątem. Chodzi o to, że oba początkowe ANDpunkty końcowe otaczającej linii muszą przecinać się z wielokątem.

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']

Wyjście:

road c
road b
road e
road f