Mehr

Finden Sie den Punkt heraus, der zwischen zwei parallelen Linien liegt


Ich habe ein Problem in ArcGIS. Ich arbeite an einer Navigationsdatenbank. In unserer Datenbank werden einspurige Straßen durch eine einzelne Linie dargestellt, während eine mehrspurige Straße (Straße mit Teiler in der Mitte) durch zwei parallele Linien (rote Linien im Bild) dargestellt wird.

Ich habe ein Punkt-Shapefile, bei dem einige Punkte innerhalb einer mehrspurigen Straße und einige außerhalb liegen.

Ich möchte ein ArcPy-Skript erstellen, das die Punkte findet, die in mehrspurige Straßen fallen. d.h. zwischen diesen parallelen Linien (im Bild markiert).

Ich weiß nicht wie ich das erreichen soll, kann mir jemand helfen

Ich habe einige Übungen damit gemacht und festgestellt, dass das Erstellen von Puffern auf einer Seite der Linie innerhalb eines mehrspurigen Polygons erstellt werden kann (im Bild gezeigt).

aber jetzt ist das Problem, dass das Polygon tatsächlich die Linie überquert (dh die mehrspurige Grenze überlappt). so wird es unnötige Punkte fangen. Gibt es eine Möglichkeit, dieses Polygon an der Straßenlinie auszurichten?

Hinweis: Integrieren funktioniert hier nicht, da es auch die Straßenlinie verschiebt. Ich muss nur das Polygon entlang der Straßenlinie ausrichten.


Ich würde es mit dem Arcpy-Algorithmus (sogar manuell!) versuchen.

  1. Finden Sie die richtige Breite der zweispurigen Straßen – hier müssen Sie möglicherweise Straßen mit der gleichen Breite gruppieren und das folgende Verfahren für jeden Cluster befolgen.
  2. Erstellen Sie Puffer beider Linien in beide Richtungen (rechts und links) mit dieser Breite (oder etwas weniger - um die Straßenfläche zu gewährleisten).
  3. Führen Sie das Schnittwerkzeug aus, um den überlappten Bereich zu erhalten.
  4. Führen Sie Nach Position auswählen aus, um Punkte auszuwählen, die innerhalb dieses Polygons liegen.

Ich würde sagen, das ist eine geometrische Übung.

PSEUDO-CODE:

  • Finden Sie für jeden Punkt (schwarzer Punkt) die nächste Straße und die Projektion des Punktes auf dieser Straße (roter Punkt).
  • Zeichnen Sie eine kurze Linie (gestrichelt) in entgegengesetzter Richtung, beginnend am schwarzen Punkt
  • Finden Sie heraus, ob es eine Kreuzung zwischen der kurzen Linie und der gleichnamigen Straße gibt, blauer Stern. Wenn es einen gibt, ist der schwarze Punkt derjenige, den wir suchen.

Wie man sieht gibt es Sonderfälle - eingekreiste schwarze Punkte:

  1. Sehr kurvenreiche 1-Linien-Straße. Dies kann vermieden werden, indem a) nur mit zweizeiligen Straßen gearbeitet wird oder b) sichergestellt wird, dass die FIDs von Straßen, die roten Punkt und Stern schneiden, unterschiedlich sind. Wenn eine kurvenreiche Straße jedoch eine Kreuzung mit einer anderen 1-Linien-Straße hat, funktioniert dies möglicherweise nicht.
  2. Der schwarze Punkt liegt auf der Verlängerung einer genau senkrechten 1-Linien-Straße. In diesem Fall besteht die Möglichkeit, dass eine einspurige Straße als nächster Nachbar ausgewählt werden kann.
  3. Schwarzer Punkt liegt auf der Linie.

Alle oben genannten Fälle sind sehr unwahrscheinlich, dennoch scheint es die sicherste Option zu sein, nur mit Straßen mit zwei Linien zu arbeiten, d. h. sie in eine separate Feature-Class zu exportieren. Fall 3 ist ein lustiger Fall, wir überlassen es dem Zufall, da die kürzeste Entfernung zur Linie nie die wahre Null ist, so dass eine 'entgegengesetzte' Richtung des Strahls gefunden werden kann, der 2 Punkte verbindet.

Python-Implementierung:

import arcpy, traceback, os, sys from arcpy import env env.overwriteoutput=True # Dinge zu ändern --------- maxD=30 mxd = arcpy.mapping.MapDocument("CURRENT") pointLR = arcpy.mapping .ListLayers(mxd,"NODES")[0] lineLR = arcpy.mapping.ListLayers(mxd,"LINKS")[0] sjOneToMany=r'D:scratchsj2.shp' RDNAME="street" # -- ----------------------- dDest=arcpy.Describe(lineLR) SR=dDest.spatialReference try: def showPyMessage(): arcpy.AddMessage(str(time .ctime()) + " - " + Nachricht) g = arcpy.Geometry() GeometrieList=arcpy.CopyFeatures_management(pointLR,g) n=len(geometryList) endPoint=arcpy.Point() arcpy.SpatialJoin_analysis(pointLR, lineLR, sjOneToMany,"JOIN_ONE_TO_MANY","KEEP_COMMON","","WITHIN_A_DISTANCE",maxD) initFidList=(-1,) für fid im Bereich(n):TARGET_FID" = %s" %str(fid) nearTable=arcpy.da .TableToNumPyArray(sjOneToMany,("TARGET_FID","JOIN_FID"),query) if len(nearTable)<2:continue fidLines=[int(row[1]) für Zeile in nearTable]FID" in %s" %str( tuple(fidLines)) listOfLines={} blackPoint=geometryLis t[fid] mit arcpy.da.SearchCursor(lineLR,("FID", "[email protected]","STREET"),query) als Zeilen: dMin=100000 für Zeile in Zeilen: shp=row[1];dCur= blackPoint.distanceTo(shp) listOfLines[row[0]]=row[-2:] if dCur

Es gibt eine andere mögliche Lösung, die vielleicht eleganter ist. Es beinhaltet Triangulation. Lassen Sie es mich wissen, wenn es von Interesse ist, und ich werde meine Antwort aktualisieren


Da die Straßen parallel verlaufen, bin ich davon ausgegangen, dass sie mit dem erstellt wurdenParallel kopierenWerkzeug in der Werkzeugleiste Bearbeiten, wodurch das Linienpaar dieselbe Richtung hat. Wir können dann über die Koordinaten der ersten Linie iterieren und sie zu einem Polygon hinzufügen und dann über die iterieren umkehren der zweiten Zeile. Es gibt definitiv einen besseren Weg, um sich an das Greifen von Linienpaaren heranzutasten; der OID-Ansatz funktioniert, aber er ist nicht sehr schön.

Sammlungen importieren arcpy FC = "fc" Punkte = "Punkte" pgons = "pgons" arcpy.env.overwriteOutput = True def buildpoly(oid_coords): #erzeuge ddict der Form OID: ddict = collections.defaultdict(list) for k,v in oid_coords: ddict[k].append(v) line1,line2 = ddict.keys() #Angenommen, die parallelen Linien haben die gleiche Richtung, also kehren Sie die zweite arr = arcpy . um .Array() arr.extend(arcpy.Point(*pt) für pt in ddict[line1]) arr.extend(arcpy.Point(*pt) für pt in ddict[line2][::-1]) return arcpy .Polygon(arr) #id ist ein ganzzahliges Feld, das parallele Linien paart unique = list(set(t[0] for t in arcpy.da.SearchCursor(FC, "id"))) polygons = [] for uni in unique: polygons.append(buildpoly([r for r in row] for row in arcpy.da.SearchCursor(FC, ["[email protected]", "[email protected]"], "id={}".format(uni) , explosion_to_points=True))) arcpy.CopyFeatures_management(polygons, pgons)

Von dort aus ist es ein Aufruf an Intersect/Select Layer nach Standort/was haben Sie. Notiere dass derSgeformtes Polygon ist nicht perfekt, da ich es freihändig gezeichnet habe und es einige Bögen gibt, dieexplodieren_to_pointsgeht nicht richtig um. Renn einfach Verdichten oder gleichwertig.