Creating a Line at a Right Angle to another line?

Question

Hello,

I have a shapefile which contains many lines features.  I would like to create a python script that will create a line and the start and end of each line, at a right angle, and a set length.

I know it will require a little trig calculations which I am happy with, but I would need to get information about the line, i.e. the cooridinates of the start and end point, the angle of the line from those points (I'm not sure how information about line shapefiles are stored)....

And then from that information create a line using python script.

Could someone point me in the right direction on how I would do this?

Thanks in advance.

FrederiqueLemieux · Answer

Thanks! I'll give this a go this week.

Will let you know how I get on.

Thanks

KimOllivier · Answer

Here is a similar problem I had a long time ago. I needed to draw a right angled line across the road at each drain (cunvert) point that was recorded in an asset management system. ArcIMS could not handle linear referencing. I got the angles by using MakeRouteEventLayer, then it was just a matter of creating a new line.

It could be a lot more elegant at 10.x but the maths is the same. I wrote out to a file and used the sample/CreateFeaturesFromTextFile, now deprecated and removed.
# culvert.py
# create culvert lines from point events
# must have loc_angle and length_m fields
# default length_m is 10 metres
# angle is from event wizard which is parallel to road
# which must be wrong at 9.1
# Input
# arg1 point event layer
# arg2 loc_angle field
# arg3 length_m field
# arg4 output featureclass
# projection from event layer
# Kim Ollivier
# 26 September 2006

import sys, string, os, win32com.client, math

def newculvert(strID,strX,strY,strAngle,strSize) :
    """
    Create a line array at right angles to angle
    and add to new featureclass already defined
    """
    id = int(strID)
    x0 = float(strX)
    y0 = float(strY)
    radAngle = float(strAngle) / 180.0 * 3.14159
    size = float(strSize) / 2.0
    if size == 0 :
        size = 10
    x1 = x0 - size * math.cos(radAngle)
    y1 = y0 - size * math.sin(radAngle)
    x2 = x0 + size * math.cos(radAngle)
    y2 = y0 + size * math.sin(radAngle)
    
    
    # print "%8d %12.3f %12.3f %8.3f %4.1f" % (id,x0,y0,angle,size)
    # print "            %12.3f %12.3f %12.3f %12.3f" % (dx1,dy1,dx2,dy2)
    print >> f1,"%d 0" % id
    print >> f1,"0 %f %f 0.0 0.0" % (x1,y1)
    print >> f1,"1 %f %f 0.0 0.0" % (x2,y2)
    return True
#----------------- main -----------------

gp = win32com.client.Dispatch("esriGeoprocessing.GpDispatch.1")

if len(sys.argv) != 6 :
    print "Need 6 parameters, setting defaults for debugging"
    strInputFeatureClass = "e:/council/clutha/shape/culvert3.shp"
    strIDField = "DRAIN_ID"
    strAngleField = "LOC_ANGLE"
    strSizeField = "LENGTH_M"
    strOutputFeatureClass = "e:/council/clutha/shape/culvert_ln.shp"
else :
    strInputFeatureClass = sys.argv[1]
    strIDField = sys.argv[2]
    strAngleField = sys.argv[3]
    strSizeField = sys.argv[4]
    strOutputFeatureClass = sys.argv[5]

if not os.path.exists(strInputFeatureClass) :
    print strInputFeatureClass,"not found"
    sys.exit()

# set the directory and workspace to the source shapefile location

os.chdir(os.path.dirname(strInputFeatureClass))
gp.Workspace = os.path.dirname(strInputFeatureClass)
f1 = open("culvert.txt","w")
print >> f1,"Polyline"

# Find the name of the shape field using Describe.
strShapeFieldName = gp.Describe(strInputFeatureClass).ShapeFieldName
# Check other fields are correct
eFlds = gp.ListFields(strInputFeatureClass)
fld = eFlds.Next()
fldOk = False
print gp.GetCount_management(strInputFeatureClass),"gp count"

while fld :
    if fld.Name.upper() == strIDField.upper() :
        fldIDOk = True
    if fld.Name.upper() == strAngleField.upper() :
        fldAngleOk = True
    if fld.Name.upper() == strSizeField.upper() :
        fldSizeOk = True
    fld = eFlds.Next()
    
if not fldIDOk:
    strErr = "No ID field found "+strIDField
    gp.AddError(strErr)
    print strErr
else :
    print strIDField," found"
if not fldAngleOk:
    strErr = "No Angle field found "+strAngleField
    gp.AddError(strErr)
    print strErr
else :
    print strAngleField," found"
if not fldSizeOk:
    strErr = "No Size field found "+strSizeField
    gp.AddError(strErr)
    print strErr
else :
    print strSizeField," found"

if not (fldIDOk or fldAngleOk or fldSizeOk) :
    print "Fatal error in field names"
    sys.exit(2)
    
nCulvert = 0
# open a Search Cursor
objRowList = gp.SearchCursor(strInputFeatureClass)
objRow = objRowList.Next()
# 
try :
    while objRow:
        objGeometry = objRow.GetValue(strShapeFieldName)
        objPart = objGeometry.GetPart(0)
        strX = objPart.X
        strY = objPart.Y
        strID = objRow.getValue(strIDField)
        strAngle = objRow.GetValue(strAngleField)
        strSize = objRow.GetValue(strSizeField)
        
        # create a new line feature
        ok = newculvert(strID,strX,strY,strAngle,strSize)
        if ok :
            nCulvert += 1
        # cleanup or looks like a memory leak
        del objRow, objPart, objGeometry
        objRow = objRowList.Next()
except Exception, e :
    print e
    gp.AddError(e)
    del objRow, objPart, objGeometry
    print "Error in Search Cursor"
    print gp.GetMessages()
    gp.AddError("Error in Search Cursor")
    sys.exit(2)

print >> f1,"END"
f1.close()

print nCulvert,"culverts total"
gp.AddMessage(str(nCulvert)+" culverts found")

# create a shapefile
gp.OverwriteOutput = 1
test = "culvert_ln.shp"
gp.CreateFeaturesFromTextFile("culvert.txt",".",test,"")
print "Created new shapefile",test
#==================================================

curtvprice · Answer

math (sorry) is a module that is part of the standard library, ie it's in all python distributions.>>> import math
>>> help(math)

If you're new to Python, be sure and check out:
Brief" rel="nofollow" target="_blank">http://docs.python.org/tutorial/stdlib.html]Brief Tour of the Standard Library