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.

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

Pingback: OSGP: Create Points at Specified Distance Interval Along a Line | Geospatiality

Hello Glen!!! Sexy blog 🙂 x

LikeLike

Hello Glen! I am an Information Systems student from Southern Brazil, I have started recently to study GIS for coastal areas voluntarily. I have been trying to write perpendicular lines all over a shapefile that contains the coastal line of Brazil, so your code has been a cricual guideline for this task which has been really difficult to find some helpful content for that task.

This code is really useful though I keep on getting an error message on line 60 when creating a new line layer.

out_ln_lyr = ds.CreateLayer(output_lns, lyr.GetSpatialRef(), ogr.wkbLineString)

it says the variable ‘lyr’ has not been declared.

I have tried to declare this variable out of the if statement to see if would solve the problem

lyr = ds.GetLayerByName(input_lyr_name)

but it didn’t work

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

I know you have written this code back in 2017 and it may be a pain in the check to check it again after a while. However, just to let you know, this piece of code is really helpful, even two years past there valuable snippets not so commonly found in python GIS tutorials. I am glad to have found your blog and I hope to be sharing my own code someday in the future.

Best Regards,

Cassiano S. Simas

LikeLike