8

I am facing one problem in ArcGIS. I work on a navigation database. In our database single lane streets are represented by a single Line, while a multi-lane street (street with divider in centre) is represented by two parallel lines (red coloured lines in picture).

I have a point shapefile with some points falling inside multi-Lane street and some outside.

I want to create an ArcPy script which would find the points which fall inside Multi-Lane Streets. i.e. between these parallel lines (marked in picture).

I do not know how to achieve this, Can somebody help me?

A multi-lane Street Problem

I did some exercise on it and i found that creating buffer on one side of line can create inside Multi-Lane polygon(Shown in Picture).

enter image description here

but now the problem is , polygon is actually crossing the line(i.e. overlapping the multi-lane boundary). so it will catch unnecessary points. is there any way to align this polygon to the street line ?

Note: integrate will not work here, because it also moves the street line. i need to just align polygon along to the street line.

  • Something like Measure the azimuth of the street - Create linestrings from each point towards angle Azimuth+90 degrees - Count how many of your parallel lines this line intersects. If zero or two -> outside, if one -> You found it. Just thinking, may work or not. Another idea is to convert the dual way street into polygon and select points which intersects is. The latter may be tricky to do with python. Well, the first as well if streets are curved. But with single sided buffer you might be able to build quite nice street polygons. – user30184 Oct 3 '15 at 18:37
  • 1
    do you have an advanced licence ? It would be quite straightforward with the near tool. – radouxju Oct 3 '15 at 18:55
  • yes i have advanced licence. – Akhil Kumar Oct 3 '15 at 19:22
  • At first i thought of taking buffer polygon and than intersecting those polygon. and find out which points fall in that intersected polygon. but the biggest problem is that In between distance is not Consistent everywhere in street. somewhere it is only 10 meter somewhere around 20 meter, in that case polygon intersect logic will be failed – Akhil Kumar Oct 3 '15 at 19:24
  • 1
    Make right side buffer of 10 m from the left side way and left side buffer from the right side one. That way you cover range 10-20 m. Overlaps do not make any harm and you can also merge the polygons first. Or make even wider one side buffer polygon and trim it by intersecting with the other side way. Use imagination and play. – user30184 Oct 3 '15 at 20:55
4

I would try below arcpy(even manual!) algorithm-

  1. Find proper width of the two lane streets- here you may need to cluster streets with same width and follow below procedure for each cluster.
  2. Create buffer both line towards both direction (right and left) with that width(or a bit less than that- to ensure road area).
  3. Run Intersection tool to get Overlapped region.
  4. Run Select by location to select points that fall inside of this polygon.
  • I think this is the way to go. Find an easy way to join the linework together, either by buffer or somehow close the lines to make a single polygon and then select within. – Barrett Oct 5 '15 at 20:35
2

I’d say this is geometric exercise.

PSEUDO CODE:

  • For every point (black point) find nearest road and find point’s projection on this road (red point).
  • Draw short line (dashed) in opposite direction starting at black point
  • Find if there is intersection between short line and same name road, blue star. If there is one, black point is the one we are after.

enter image description here

As one can see there are special cases – circled black points:

  1. Very twisty 1 line road. This can be eliminated by a) working with 2 line roads only or b) making sure FIDs of roads that intersect red dot and star are different. However if bendy road has a junction with another 1 line road, this might not work.
  2. Black point is sitting on extension of exactly perpendicular 1 line road. In this case there is a chance that 1 lane road can be picked as a nearest neighbour.
  3. Black point sits on the line.

All of above cases are very unlikely, nevertheless it seems that the safest option is to work with 2 line roads only, i.e. export them to a separate feature class. Case 3 is a funny one, we’ll leave it to chance, because shortest distance to line is never true zero, thus ‘opposite’ direction of ray connecting 2 points can be found.

Python implementation:

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()            

There is another possible solution perhaps more elegant. It involves triangulation. Let me know if it is of interest and I'll update my answer

  • This is quite complex, wow. It seems like it would be a lot simpler to create a polygon from the lines and then use ray-casting. Determining if a point is on a line should be straightforward as well. – Paul Oct 5 '15 at 21:10
  • 1
    If you are able to create polygons from correct lines, no need in casting. Select by location will do. Creating polygons is a challenge though – FelixIP Oct 5 '15 at 21:12
  • Will it work well at bends- just for clarification:) – SIslam Oct 6 '15 at 5:07
  • 1
    @SIslam it should work even with big bends similar to case 1 (see if n==fidNear:continue) line. Well, if there is no 1 lane road going in. I keep thinking that dissolve can help, but not always – FelixIP Oct 6 '15 at 5:22
  • @Islam Oops! It will not, because condition (if n==fidNear:continue)eliminates points sitting outside the bend, but mark point inside as one sitting outside. Sharp turn needed though, radius smaller than width? – FelixIP Oct 6 '15 at 5:36
0

Since the streets are parallel, I've assumed that they were created with the Copy Parallel tool in the Edit toolbar thus making the pair of lines have the same direction. We can then iterate over the coordinates of the first line and add them to a polygon and then iterate over the reverse of the second line. There's definitely a better way to approach grabbing line pairs; the OID approach works, but it isn't very pretty.

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)

From there, it's a call to Intersect/Select Layer by location/what have you. Note that the S shaped polygon isn't perfect since I drew it freehand and there are some arcs which explode_to_points doesn't handle properly. Just run Densify or equivalent.

enter image description here

  • This is road network dataset, thus 1 lane roads are connected to 2 lanes via node, i.e. there is no such things as pairs of parallel features – FelixIP Oct 5 '15 at 22:47
  • You might want to extend your solution adding dissolving by individual road names first (no m-parts) and consider cases of 1 or 2 lines as a result – FelixIP Oct 5 '15 at 22:57
  • @FelixIP, I'm not very familiar with network datasets. My solution was mostly a proof of concept of how it can be done with simple lines (OP can extend it to cover m resolution, multipart, etc.). I don't know how features like this are actually represented in a network. – Paul Oct 5 '15 at 23:04
  • @ Paul Same name road can be represented by 100s segments sitting in different rows in the table. Moreover double lane road could become single lane somewhere. Dissolve will fail badly if no of parts not in (1,2), this is why I didn’t go with triangulation solution – FelixIP Oct 5 '15 at 23:12
  • 1
    @AkhilKumar, it does not matter if they are roughly parallel. This is tracing the existing lines. – Paul Oct 6 '15 at 16:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.