Trouver le point qui se situe entre deux lignes parallèles

Je rencontre un problème dans ArcGIS. Je travaille sur une base de données de navigation. Dans notre base de données, les rues à voie unique sont représentées par une seule ligne, tandis qu'une rue à plusieurs voies (rue avec séparateur au centre) est représentée par deux lignes parallèles (lignes rouges sur l'image).

J'ai un fichier de formes de points avec certains points tombant à l'intérieur d'une rue à plusieurs voies et d'autres à l'extérieur.

Je veux créer un script ArcPy qui trouverait les points qui se trouvent dans les rues à plusieurs voies. c'est-à-dire entre ces lignes parallèles (marquées sur l'image).

Je ne sais pas comment y parvenir, quelqu'un peut-il m'aider?

J'ai fait un peu d'exercice dessus et j'ai trouvé que la création d'un tampon sur un côté de la ligne peut créer à l'intérieur du polygone à plusieurs voies (montré sur l'image).

mais maintenant le problème est que le polygone traverse réellement la ligne (c'est-à-dire qu'il chevauche la frontière à plusieurs voies). il capturera donc des points inutiles. existe-t-il un moyen d'aligner ce polygone sur la ligne de rue?

Remarque: l'intégration ne fonctionnera pas ici, car elle déplace également la ligne de rue. j'ai juste besoin d'aligner le polygone le long de la ligne de rue.

J'essaierais ci-dessous l'algorithme arcpy (même manuel!)

1. Trouvez la bonne largeur des rues à deux voies - ici, vous devrez peut-être regrouper des rues de même largeur et suivre la procédure ci-dessous pour chaque groupe.
2. Créez un tampon les deux lignes dans les deux sens (droite et gauche) avec cette largeur (ou un peu moins que cela - pour assurer la zone de la route).
3. Exécutez l'outil Intersection pour obtenir la région chevauchée.
4. Exécutez Sélectionner par emplacement pour sélectionner les points qui se trouvent à l'intérieur de ce polygone.

Je dirais que c'est un exercice géométrique.

CODE PSEUDO:

• Pour chaque point (point noir), trouvez la route la plus proche et trouvez la projection du point sur cette route (point rouge).
• Tracez une ligne courte (en pointillés) dans la direction opposée à partir du point noir
• Trouvez s'il y a une intersection entre la ligne courte et la route du même nom, l'étoile bleue. S'il y en a un, le point noir est celui que nous recherchons.

Comme on peut le voir, il existe des cas particuliers - points noirs encerclés:

1. Route 1 ligne très sinueuse. Ceci peut être éliminé en a) travaillant uniquement avec des routes à 2 lignes ou b) en s'assurant que les FID des routes qui croisent le point rouge et l'étoile sont différents. Cependant, si la route coudée a une jonction avec une autre route à 1 ligne, cela pourrait ne pas fonctionner.
2. Le point noir est assis sur le prolongement d'une route à 1 ligne exactement perpendiculaire. Dans ce cas, il est possible que la route à 1 voie soit choisie comme voisin le plus proche.
3. Le point noir se trouve sur la ligne.

Tous les cas ci-dessus sont très improbables, mais il semble que l'option la plus sûre consiste à travailler uniquement avec des routes à 2 lignes, c'est-à-dire à les exporter vers une classe d'entités distincte. Le cas 3 est drôle, nous allons le laisser au hasard, car la distance la plus courte par rapport à la ligne n'est jamais vraiment nulle, donc la direction `` opposée '' du rayon reliant 2 points peut être trouvée.

Implémentation de Python:

import arcpy, traceback, os, sys
from arcpy import env
env.overwriteoutput=True

# things to change ---------
maxD=30
mxd = arcpy.mapping.MapDocument("CURRENT")
pointLR = arcpy.mapping.ListLayers(mxd,"NODES")[0]
lineLR = arcpy.mapping.ListLayers(mxd,"LINKS")[0]
sjOneToMany=r'D:\scratch\sj2.shp'
RDNAME='street'
# -------------------------
dDest=arcpy.Describe(lineLR)
SR=dDest.spatialReference

try:
    def showPyMessage():
        arcpy.AddMessage(str(time.ctime()) + " - " + message)
    g = arcpy.Geometry()
    geometryList=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,)
    for fid in range(n):
        query='"TARGET_FID" = %s' %str(fid)
        nearTable=arcpy.da.TableToNumPyArray(sjOneToMany,("TARGET_FID","JOIN_FID"),query)
        if len(nearTable)<2:continue
        fidLines=[int(row[1]) for row in nearTable]
        query='"FID" in %s' %str(tuple(fidLines))
        listOfLines={}
        blackPoint=geometryList[fid]
        with arcpy.da.SearchCursor(lineLR,("FID", "Shape@","STREET"),query) as rows:
            dMin=100000
            for row in rows:
                shp=row[1];dCur=blackPoint.distanceTo(shp)
                listOfLines[row[0]]=row[-2:]
                if dCur<dMin:
                    fidNear,lineNear, roadNear=row
                    dMin=dCur
            chainage=lineNear.measureOnLine(blackPoint)
            redPoint=lineNear.positionAlongLine (chainage).firstPoint
            smallD=blackPoint.distanceTo(redPoint)
            fp=blackPoint.firstPoint
            dX=(redPoint.X-fp.X)*(maxD-smallD)/smallD
            dY=(redPoint.Y-fp.Y)*(maxD-smallD)/smallD
            endPoint.X=fp.X-dX;endPoint.Y=fp.Y-dY
            dashLine=arcpy.Polyline(arcpy.Array([fp,endPoint]),SR)

            for n in listOfLines:
                if n==fidNear:continue
                line, road=listOfLines[n]
                if road!=roadNear:continue
                blueStars=dashLine.intersect(line,1)
                if blueStars.partCount==0:continue
                initFidList+=(fid,); break
    query='"FID" in %s' %str(initFidList)
    arcpy.SelectLayerByAttribute_management(pointLR, "NEW_SELECTION", query)
    arcpy.AddMessage ('\n %i point(s) found' %(len(initFidList)-1))
except:
    message = "\n*** PYTHON ERRORS *** "; showPyMessage()
    message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyMessage()
    message = "Python Error Info: " +  str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()            

Il existe une autre solution possible, peut-être plus élégante. Cela implique une triangulation. Faites-moi savoir si cela vous intéresse et je mettrai à jour ma réponse

Étant donné que les rues sont parallèles, j'ai supposé qu'elles ont été créées avec l' Copy Paralleloutil dans la barre d'outils Modifier, ce qui fait que la paire de lignes a la même direction. Nous pouvons ensuite itérer sur les coordonnées de la première ligne et les ajouter à un polygone, puis itérer sur l' inverse de la deuxième ligne. Il y a certainement une meilleure façon d'approcher les paires de lignes de saisie; l'approche OID fonctionne, mais elle n'est pas très jolie.

import collections
import arcpy

FC = "fc"
points = "points"
pgons = "pgons"
arcpy.env.overwriteOutput = True

def buildpoly(oid_coords):
    #create ddict of the form OID:<x1y1, x2y2, ..., xn-1yn-1, xnyn>
    ddict = collections.defaultdict(list)    
    for k,v in oid_coords:
        ddict[k].append(v)

    line1,line2 = ddict.keys()    

    #Assume that the parallel lines have same direction, so reverse the second
    arr = arcpy.Array()
    arr.extend(arcpy.Point(*pt) for pt in ddict[line1])    
    arr.extend(arcpy.Point(*pt) for pt in ddict[line2][::-1])

    return arcpy.Polygon(arr)

#id is an integer field that pairs parallel lines together
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,
                                                                                ["OID@", "SHAPE@XY"],
                                                                                "id={}".format(uni),
                                                                                explode_to_points=True)))


arcpy.CopyFeatures_management(polygons, pgons)

À partir de là, c'est un appel à Intersect / Select Layer by location / what have you. Notez que le Spolygone formé n'est pas parfait car je l'ai dessiné à main levée et il y a des arcs qui explode_to_pointsne se gèrent pas correctement. Exécutez simplement Densifyou équivalent.