Tag Archives: Distance

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

Posted on by
3

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
    bearing = math.radians(angle)
    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):
    bearing = math.radians(bearing)
    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
        shapely_line = MultiLineString(wkt.loads(line_geom))
        ## 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.

chainage_ticks_along_linechainage_ticks_along_line_attributes

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

OSGP: Create Points at Specified Distance Interval Along a Line

Posted on by
3

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:
    print "{0} NOT found in {1}".format(input_lyr_name, gdb)
    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
        shapely_line = MultiLineString(wkt.loads(line_geom))
        ## 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)
            pnt.AddPoint(pt.x, pt.y)
            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.

distance_points_along_linedistance_points_along_line_attributesPlease 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.