# OSGP: Create Chainage ticks along a Line at Specified Distance Intervals

This builds on from the previous post creating points at specified distances along a line. Here, we create perpendicular chainage ticks that traverse the main line.

```from osgeo import ogr
from shapely.geometry import MultiLineString, LineString, Point
from shapely import wkt
import sys, math

## http://wikicode.wikidot.com/get-angle-of-line-between-two-points
## angle between two points
def getAngle(pt1, pt2):
x_diff = pt2.x - pt1.x
y_diff = pt2.y - pt1.y
return math.degrees(math.atan2(y_diff, x_diff))

## start and end points of chainage tick
## get the first end point of a tick
def getPoint1(pt, bearing, dist):
angle = bearing + 90
x = pt.x + dist * math.cos(bearing)
y = pt.y + dist * math.sin(bearing)
return Point(x, y)
## get the second end point of a tick
def getPoint2(pt, bearing, dist):
x = pt.x + dist * math.cos(bearing)
y = pt.y + dist * math.sin(bearing)
return Point(x, y)

## set the driver for the data
driver = ogr.GetDriverByName("FileGDB")
## path to the FileGDB
gdb = r"C:\Users\******\Documents\ArcGIS\Default.gdb"
## open the GDB in write mode (1)
ds = driver.Open(gdb, 1)

## linear feature class
input_lyr_name = "input_line"

## distance between each points
distance = 10
## the length of each tick
tick_length = 20

## output tick line fc name
output_lns = "{0}_{1}m_lines".format(input_lyr_name, distance)

## list to hold all the point coords
list_points = []

## reference the layer using the layers name
if input_lyr_name in [ds.GetLayerByIndex(lyr_name).GetName() for lyr_name in range(ds.GetLayerCount())]:
lyr = ds.GetLayerByName(input_lyr_name)
print "{0} found in {1}".format(input_lyr_name, gdb)

## if the output already exists then delete it
if output_lns in [ds.GetLayerByIndex(lyr_name).GetName() for lyr_name in range(ds.GetLayerCount())]:
ds.DeleteLayer(output_lns)
print "Deleting: {0}".format(output_lns)

## create a new line layer with the same spatial ref as lyr
out_ln_lyr = ds.CreateLayer(output_lns, lyr.GetSpatialRef(), ogr.wkbLineString)

## distance/chainage attribute
chainage_fld = ogr.FieldDefn("CHAINAGE", ogr.OFTReal)
out_ln_lyr.CreateField(chainage_fld)
## check the geometry is a line
first_feat = lyr.GetFeature(1)

## accessing linear feature classes using FileGDB driver always returns a MultiLinestring
if first_feat.geometry().GetGeometryName() in ["LINESTRING", "MULTILINESTRING"]:
for ln in lyr:
## list to hold all the point coords
list_points = []
## set the current distance to place the point
current_dist = distance
## get the geometry of the line as wkt
line_geom = ln.geometry().ExportToWkt()
## make shapely MultiLineString object
## get the total length of the line
line_length = shapely_line.length
## append the starting coordinate to the list
list_points.append(Point(list(shapely_line[0].coords)[0]))
## https://nathanw.net/2012/08/05/generating-chainage-distance-nodes-in-qgis/
## while the current cumulative distance is less than the total length of the line
while current_dist < line_length:
## use interpolate and increase the current distance
list_points.append(shapely_line.interpolate(current_dist))
current_dist += distance
## append end coordinate to the list
list_points.append(Point(list(shapely_line[0].coords)[-1]))

## add lines to the layer
## this can probably be cleaned up better
## but it works and is fast!
for num, pt in enumerate(list_points, 1):
## start chainage 0
if num == 1:
angle = getAngle(pt, list_points[num])
line_end_1 = getPoint1(pt, angle, tick_length/2)
angle = getAngle(line_end_1, pt)
line_end_2 = getPoint2(line_end_1, angle, tick_length)
tick = LineString([(line_end_1.x, line_end_1.y), (line_end_2.x, line_end_2.y)])
feat_dfn_ln = out_ln_lyr.GetLayerDefn()
feat_ln = ogr.Feature(feat_dfn_ln)
feat_ln.SetGeometry(ogr.CreateGeometryFromWkt(tick.wkt))
feat_ln.SetField("CHAINAGE", 0)
out_ln_lyr.CreateFeature(feat_ln)

## everything in between
if num < len(list_points) - 1:
angle = getAngle(pt, list_points[num])
line_end_1 = getPoint1(list_points[num], angle, tick_length/2)
angle = getAngle(line_end_1, list_points[num])
line_end_2 = getPoint2(line_end_1, angle, tick_length)
tick = LineString([(line_end_1.x, line_end_1.y), (line_end_2.x, line_end_2.y)])
feat_dfn_ln = out_ln_lyr.GetLayerDefn()
feat_ln = ogr.Feature(feat_dfn_ln)
feat_ln.SetGeometry(ogr.CreateGeometryFromWkt(tick.wkt))
feat_ln.SetField("CHAINAGE", distance * num)
out_ln_lyr.CreateFeature(feat_ln)

## end chainage
if num == len(list_points):
angle = getAngle(list_points[num - 2], pt)
line_end_1 = getPoint1(pt, angle, tick_length/2)
angle = getAngle(line_end_1, pt)
line_end_2 = getPoint2(line_end_1, angle, tick_length)
tick = LineString([(line_end_1.x, line_end_1.y), (line_end_2.x, line_end_2.y)])
feat_dfn_ln = out_ln_lyr.GetLayerDefn()
feat_ln = ogr.Feature(feat_dfn_ln)
feat_ln.SetGeometry(ogr.CreateGeometryFromWkt(tick.wkt))
feat_ln.SetField("CHAINAGE", int(line_length))
out_ln_lyr.CreateFeature(feat_ln)

del ds```

The output is a linear feature class containing chainage ticks and distance attribute.

Please leave comments if this can be improved or if you found it useful.

# OSGP: Create Points at Specified Distance Interval Along a Line

This workflow with Python using OGR and Shapely creates points along a line at a specified distance interval. I use the FileGDB driver here to read from and write data to but you can change these to suit your requirements. The code is commented…

```from osgeo import ogr
from shapely.geometry import MultiLineString, Point
from shapely import wkt
import sys

## set the driver for the data
driver = ogr.GetDriverByName("FileGDB")
################################################################################
## CHANGE gdb, input_lyr_name, distance, output_pts (optional)

## path to the FileGDB
gdb = r"C:\Users\******\Documents\ArcGIS\Default.gdb"
## open the GDB in write mode (1)
ds = driver.Open(gdb, 1)

## single linear feature
input_lyr_name = "input_line"

## distance between each points
distance = 10

## output point fc name
output_pts = "{0}_{1}m_points".format(input_lyr_name, distance)
################################################################################

## reference the layer using the layers name
if input_lyr_name in [ds.GetLayerByIndex(lyr_name).GetName() for lyr_name in range(ds.GetLayerCount())]:
lyr = ds.GetLayerByName(input_lyr_name)
print "{0} found in {1}".format(input_lyr_name, gdb)
## if the feature class cannot be found exit gracefully
else:
sys.exit()

## if the output already exists then delete it
if output_pts in [ds.GetLayerByIndex(lyr_name).GetName() for lyr_name in range(ds.GetLayerCount())]:
ds.DeleteLayer(output_pts)
print "Deleting: {0}".format(output_pts)

## create a new point layer with the same spatial ref as lyr
out_lyr = ds.CreateLayer(output_pts, lyr.GetSpatialRef(), ogr.wkbPoint)

## create a field to hold the distance values
dist_fld = ogr.FieldDefn("DISTANCE", ogr.OFTReal)
out_lyr.CreateField(dist_fld)
## check the geometry is a line
first_feat = lyr.GetFeature(1)

## accessing linear feature classes using FileGDB driver always returns a MultiLinestring
if first_feat.geometry().GetGeometryName() in ["LINESTRING", "MULTILINESTRING"]:
for ln in lyr:
## list to hold all the point coords
list_points = []
## set the current distance to place the point
current_dist = distance
## get the geometry of the line as wkt
line_geom = ln.geometry().ExportToWkt()
## make shapely MultiLineString object
## get the total length of the line
line_length = shapely_line.length
## append the starting coordinate to the list
list_points.append(Point(list(shapely_line[0].coords)[0]))
## https://nathanw.net/2012/08/05/generating-chainage-distance-nodes-in-qgis/
## while the current cumulative distance is less than the total length of the line
while current_dist < line_length:
## use interpolate and increase the current distance
list_points.append(shapely_line.interpolate(current_dist))
current_dist += distance
## append end coordinate to the list
list_points.append(Point(list(shapely_line[0].coords)[-1]))

## add points to the layer
## for each point in the list
for num, pt in enumerate(list_points, 1):
## create a point object
pnt = ogr.Geometry(ogr.wkbPoint)
feat_dfn = out_lyr.GetLayerDefn()
feat = ogr.Feature(feat_dfn)
feat.SetGeometry(pnt)
## populate the distance values for each point.
## start point
if num == 1:
feat.SetField("DISTANCE", 0)
elif num < len(list_points):
feat.SetField("DISTANCE", distance * (num - 1))
## end point
elif num == len(list_points):
feat.SetField("DISTANCE", int(line_length))
## add the point feature to the output.
out_lyr.CreateFeature(feat)

else:
print "Error: make sure {0} is a linear feature class with at least one feature".format(input_lyr_name)
sys.exit()

del ds, out_lyr```

This will create a point feature class with a single attribute containing the distance value from the start of the line.

Please leave any constructive feedback if you think this can be improved or if it worked for you.

Also see Create Chainage Ticks for creating perpendicular lines traversing the main line at specified distances.