Creating perpendicular line from points to existing line?

Question

I would like to create perpendicular lines from points to line, similar to this:

I have ArcGIS 10.1, so "queryPointAndDistance" doesn´t work for me, and basic licence level, so I can´t use "Generate Near Table" neither.

I thought about following approach:

• Coordinates X1Y1 are known (coordinates of the point feature)
• acquire X2Y2 coordinates of the closest point on the line segment
• create line feature from X1Y1 to X2Y2
• new line feature is perpendicular to existing line

I tried to use "distanceTo" method from point feature to line feature, but I was unable to aquire coordinates of the closest point of the line segment (X2Y2), which should be used during "distanceTo" method. Is it possible to get coordinates of the closest point on the line segment from "distanceTo" method or should I use another approach?

point_feature = "points_select"
line_feature = "line_select"

pointcursor = arcpy.da.SearchCursor(point_feature, ['SHAPE@'])        

for row in pointcursor:
      geometry = pointcursor[0]
      linecursor = arcpy.da.SearchCursor(line_feature, ['SHAPE@'])
      for row in linecursor:
          newgeometry = linecursor[0]
          dist = newgeometry.distanceTo(geometry)
      print dist

Answer 1

It took below script 3 seconds to find 100 nearest points on the line. It uses positionAlongLine method on geometry which I believe available in 10.1

SCRIPT:

import  arcpy, math
from arcpy import env
env.overwriteOutput = True
# LAYERS IN TOC OF MXD
points,pipes,connections="PARC_POINTS","PIPES","CONNECTIONS"
sj=r"in_memory\sj"
# golden section method to find minimum of function below
gr,hypot=(math.sqrt(5)-1)/2,math.hypot
def gss(a,b,tol):
    c=b-gr*(b-a)
    d=a+gr*(b-a)
    while abs(c-d)>tol:       
        fc=f(c);fd=f(d)
        if fc<fd:
            b,d=d,c
            c=b-gr*(b-a)
        else:
            a,c=c,d
            d=a+gr*(b-a)
    return (b+a)/2
# function to minimise
def f(x):
    theP=pipe.positionAlongLine(x).firstPoint
    dist=hypot(x1-theP.X,y1-theP.Y)
    return dist
# find pipe nearest to point. NOTE: field PIPE_ID in pipes is INTEGER!
arcpy.SpatialJoin_analysis(points, pipes, sj, "JOIN_ONE_TO_ONE",
                           "KEEP_ALL", """PIPE_ID "PIPE_ID" true true false 5 Short 0 5 ,First,#,PIPES,PIPE_ID,-1,-1""",
                           "CLOSEST")
curT=arcpy.da.InsertCursor(connections,"Shape@")
# shuffle through points
with arcpy.da.SearchCursor(sj,["Shape@","PIPE_ID"]) as cursor:
    for pG,pipeID in cursor:
        p1=pG.firstPoint
        x1,y1=p1.X,p1.Y
        dq ='"PIPE_ID"=%s'%pipeID
        with arcpy.da.SearchCursor(pipes,"Shape@", dq) as pipeCursor:
            for row in pipeCursor:pipe=row[0];break
        L=pipe.length
# find nearest point on pipe (accurate to 1 cm, see last value)
        chainage=gss(0.0,L,0.01)
        p2=pipe.positionAlongLine(chainage).firstPoint
# and connect original point to it by line
        arr=arcpy.Array([p1,p2])
        perpendicular=arcpy.Polyline(arr)
        curT.insertRow((perpendicular,))

OUTPUT:

I hope script is documented well enough to understand naming conventions and must have field.

Answer 2

I tryed once more (sorry about non english notes):

import arcpy
arcpy.env.overwriteOutput = True

#input:
line_upper = arcpy.GetParameterAsText(0) #input shape # Shapefile
line_lower = arcpy.GetParameterAsText(1) #input shape # Shapefile

step = arcpy.GetParameter(2) #step #double
step = int(step)

max_distance_plus = arcpy.GetParameter(3) #step #double
max_distance_plus = int(max_distance_plus)

#input edited
max_distance_minus = max_distance_plus *(-1)

desc_linie_upper = arcpy.Describe(line_upper) 
line_upper_path = desc_linie_upper.path + "\\"

## temp
pracovni_path = line_upper_path + "__pracovni_smazat" + "\\"

if not arcpy.Exists(pracovni_path):
      arcpy.CreateFolder_management(line_upper_path, "__pracovni_smazat")

Route_ID = "ID"
txt_plus_path = pracovni_path + "_test_plus.txt"
txt_minus_path = pracovni_path + "_test_minus.txt"

PointsAlongLines_fc = "_PointsAlongLines.shp"
PointsAlongLines = pracovni_path + "_PointsAlongLines.shp"
RouteEvent_merge = pracovni_path + "_RouteEvent_merge.shp"
PointsToLine = pracovni_path + "_PointsToLine.shp"
IntersectPoints = pracovni_path + "_IntersectPoints.shp"
PointsMerge = pracovni_path + "_pointsMerge.shp"
IntersectPoints_single_fc = "_IntersectPoints_single.shp"
IntersectPoints_single = pracovni_path + "_IntersectPoints_single.shp"
MergeIntersectPpal = pracovni_path + "_MergeIntersectPpal.shp"

#fin output file
svah_feature_fc = "svahove_linie.shp"
svah_feature = line_upper_path + svah_feature_fc

#txt
PlusWrite = open(txt_plus_path,"w")
hlavicka_sloupcu = "Location" + "\t" + "Offset" +"\t"+ "Route_ID"
PlusWrite.writelines(hlavicka_sloupcu + "\n")

MinusWrite = open(txt_minus_path,"w")
hlavicka_sloupcu = "Location" + "\t" + "Offset" +"\t"+ "Route_ID"
MinusWrite.writelines(hlavicka_sloupcu + "\n")

#Position along lines
ppal_all = []
with arcpy.da.SearchCursor(line_upper, ["SHAPE@",Route_ID]) as ppcur:
      for pprow in ppcur:
            length = int(pprow[0].length)
            for i in xrange(0, length + step, step):
                  ppal = pprow[0].positionAlongLine(i)
                  ppal_all.append([ppal,i])

                  #zapis do textaku
                  PlusWrite.writelines(str(i) + "\t" + str(max_distance_plus) + "\t" + str(pprow[1]) +"\n")
                  MinusWrite.writelines(str(i) + "\t" + str(max_distance_minus) + "\t" + str(pprow[1]) +"\n")
del ppcur
PlusWrite.close()
MinusWrite.close()

#points along line new shp
arcpy.CreateFeatureclass_management(pracovni_path,PointsAlongLines_fc,"POINT")
arcpy.AddField_management(PointsAlongLines,"Location","DOUBLE")

with arcpy.da.InsertCursor(PointsAlongLines, ["SHAPE@","Location"]) as icur:
      for i in ppal_all:
            icur.insertRow(i)           
del icur        

#create line from points along line buffer
arcpy.MakeRouteEventLayer_lr(line_upper, Route_ID, txt_plus_path,"Route_ID POINT Location","RouteEvent_Plus", "Offset")
arcpy.MakeRouteEventLayer_lr(line_upper, Route_ID, txt_minus_path,"Route_ID POINT Location","RouteEvent_Minus", "Offset")
arcpy.Merge_management(["RouteEvent_Minus","RouteEvent_Plus"],RouteEvent_merge)

arcpy.PointsToLine_management(RouteEvent_merge,PointsToLine,"Location")

#Intersect
arcpy.Intersect_analysis([PointsToLine,line_lower],IntersectPoints,"ALL",0.01,"POINT")

IntersectPoints_list = []
with arcpy.da.SearchCursor(IntersectPoints, ["SHAPE@","Location"]) as iscur:
      for isrow in iscur:            
            for part in isrow[0]:                  
                  IntersectPoints_list.append([part,part.M])

del iscur

#multi to single points
arcpy.CreateFeatureclass_management(pracovni_path,IntersectPoints_single_fc,"POINT")
arcpy.AddField_management(IntersectPoints_single,"Location","DOUBLE")

with arcpy.da.InsertCursor(IntersectPoints_single, ["SHAPE@","Location"]) as icur:
      for i in IntersectPoints_list:
            icur.insertRow(i)           
del icur      

arcpy.CreateFeatureclass_management(line_upper_path,IntersectPoints_single_fc,"POINT")
arcpy.AddField_management(IntersectPoints_single,"Location","DOUBLE")

#fin line
exp = arcpy.AddFieldDelimiters(IntersectPoints_single,"Location")
InterPal_list =[]
with arcpy.da.SearchCursor(IntersectPoints_single,["SHAPE@","Location"]) as iscur:
      for isrow in iscur:
            isrowgeom = isrow[0]            
            with arcpy.da.SearchCursor(PointsAlongLines,["SHAPE@","Location","FID"], exp + " = " + str(isrow[1])) as pscur:
                  for psrow in pscur:                        
                        psrowgeom = psrow[0]
                        if psrow[2] % 2 ==0:                              
                              InterPal_list.append(arcpy.Polyline(arcpy.Array([isrowgeom.firstPoint,psrowgeom.firstPoint])))
                        else:
                              npnt = arcpy.Point()
                              npnt.X = ((isrowgeom.firstPoint.X - psrowgeom.firstPoint.X)/4+psrowgeom.firstPoint.X)
                              npnt.Y = ((isrowgeom.firstPoint.Y - psrowgeom.firstPoint.Y)/4+psrowgeom.firstPoint.Y)
                              InterPal_list.append(arcpy.Polyline(arcpy.Array([npnt,psrowgeom.firstPoint])))
del iscur
del pscur

##new/existing shape
if arcpy.Exists(svah_feature):
      arcpy.Append_management(InterPal_list,svah_feature,"NO_TEST")
      arcpy.AddMessage("\n" + "shape: " + svah_feature + "\n")
else:
      arcpy.CopyFeatures_management(InterPal_list, svah_feature)
      arcpy.AddMessage("\n" + "novy shape: " + svah_feature + "\n")

Answer 3

I mentioned before, creating perpendicular lines is part of the bigger question. I tried to create slope hatches. I think I partly succeeded, It works only for pair of lines, upper and lower. Surely there is easier and more accurate way, but for now it's enough.

SCRIPT:

import arcpy 

#parameters
line_upper = arcpy.GetParameterAsText(0) #input shape # Shapefile
line_lower = arcpy.GetParameterAsText(1) #input shape # Shapefile
step = arcpy.GetParameter(2) #step #double
step = int(step)

#Position along lines
PointsAlongLines = "in_memory/PointsAlongLines"

ppal_all = []
with arcpy.da.SearchCursor(line_upper, ["SHAPE@"]) as ppcur:
      for pprow in ppcur:
            length = int(pprow[0].length)
            for i in xrange(0, length + step, step):
                  ppal = pprow[0].positionAlongLine(i)
                  ppal_all.append(ppal)
arcpy.CopyFeatures_management(ppal_all,PointsAlongLines)

#clean
del ppcur

#parameters
svah_feature = "svahove_linie.shp"
desc_linie_upper = arcpy.Describe(line_upper) 
line_upper_path = desc_linie_upper.path + "\\"
buff_feature_path = line_upper_path + "_buff_feature.shp"
svah_feature_path = line_upper_path + svah_feature

poradi = 0
#min distance line - point
arcpy.AddField_management(PointsAlongLines,"Min_dist","DOUBLE")
with arcpy.da.UpdateCursor(PointsAlongLines, ["SHAPE@", "Min_dist"]) as pcur:
      for prow in pcur:
            point_geometry = pcur[0]            
            with arcpy.da.SearchCursor(line_lower, ["SHAPE@"], ) as lcur:
                  for lrow in lcur:
                        line_geometry = lcur[0]
                        mindist = line_geometry.distanceTo(point_geometry)
                        poradi = poradi+1
                        if poradi % 2 ==0:
                              prow[1] = mindist/2
                        else:
                              prow[1] = mindist                                                                                   
                        pcur.updateRow(prow)
#clean                        
del pcur
del lcur

#buffer 
arcpy.Buffer_analysis(PointsAlongLines, buff_feature_path,"Min_dist")

svah_ary = []
#vertex by vertex
with arcpy.da.SearchCursor(buff_feature_path, ["SHAPE@","SHAPE@TRUECENTROID"]) as bcur:
      for brow in bcur:            
            buff_list = []
            for part in brow[0]:                  
                  for pnt in part:
                        #list min distances
                        with arcpy.da.SearchCursor(line_lower, ["SHAPE@"], ) as lcur:
                              for lrow in lcur:
                                    line_geometry = lcur[0]
                                    buff_line_dist = line_geometry.distanceTo(arcpy.Point(pnt.X,pnt.Y))                                    
                                    buff_list.append([pnt.X, pnt.Y,buff_line_dist])
                  #sorting list
                  buff_list = sorted(buff_list, key=lambda buff_list_entry: buff_list_entry[2])
                  #array
                  svah_ary.append(arcpy.Array([arcpy.Point(brow[1][0], brow[1][1]), arcpy.Point(buff_list[0][0], buff_list[0][1])]))                          

#clean                  
del lcur
del bcur

#new/existing shape
svahove_linie = []
for part in svah_ary:
    svahove_linie.append(arcpy.Polyline(arcpy.Array(part)))
if arcpy.Exists(svah_feature_path):
      arcpy.Append_management(svahove_linie,svah_feature_path,"NO_TEST")
      arcpy.AddMessage("\n" + "shape: " + svah_feature_path + "\n")
else:
      arcpy.CopyFeatures_management(svahove_linie, svah_feature_path)
      arcpy.AddMessage("\n" + "novy shape: " + svah_feature_path + "\n")

#clean   
arcpy.Delete_management(buff_feature_path)            
arcpy.Delete_management(PointsAlongLines)