xander_bakker

Translating the lowest overhead power lines into a continuous surface

Blog Post created by xander_bakker on Jul 18, 2016

Let’s try and tackle the problem of translating the lowest overhead power lines into a continuous surface. We could use the lines and create a TIN, but we want to have some information outside the power lines, since there may be trees there that interfere with the power line. If you edit the featureclass and create a line, using an offset, the new geometry will be 2D, since the Z value will not be passed to the new geometry. So that’s not an option. Probably the easiest way is to achieve this is using some python code.

 

What we want to achieve is creating points perpendicular to the power lines, with the Z value of the overhead power line. Although we could use the method described in my blog What if you need to create a 3D line representing the movement due to wind of an overhead power lines?  to generate points or lines and generate a TIN from them.

 

Code used to generate perpendicular points to 3D line:

 

import arcpy

def main():

    # configuracion
    fc = r'D:\Xander\_LiDAR\ElSalto\tmp.gdb\vanos'
    fc_out = r'D:\Xander\_LiDAR\ElSalto\tmp.gdb\cable_pnts_Z_v01'
    sr = arcpy.Describe(fc).spatialReference
    d_max = 40
    step_size = 5
    where = 'vano = 1'

    # cursor para leer lineas 3D de cables
    lst_pntg = []
    with arcpy.da.SearchCursor(fc, ('SHAPE@'), where) as curs:
        # para cada feature
        for row in curs:
            line = row[0]
            # definir angulo del vano
            angle = getAngle(line, sr)

            #  para cada parte de la linea
            for part in line:
                # para cada punto en la parte actual
                for pnt in part:
                    pntg = arcpy.PointGeometry(pnt, sr)
                    z = pnt.Z
                    # crear puntos perpendiculares
                    for i in range(int(2 * d_max/step_size)+1):
                        d = i * step_size - d_max
                        pntg_tmp1 = pntg.pointFromAngleAndDistance(angle-90, d, 'PLANAR')
                        pntg_tmp2 = pntg.pointFromAngleAndDistance(angle+90, d, 'PLANAR')
                        pnt1 = arcpy.Point(pntg_tmp1.firstPoint.X, pntg_tmp1.firstPoint.Y, z)
                        pnt2 = arcpy.Point(pntg_tmp2.firstPoint.X, pntg_tmp2.firstPoint.Y, z)
                        pntg_out1 = arcpy.PointGeometry(pnt1, sr, True)
                        pntg_out2 = arcpy.PointGeometry(pnt2, sr, True)
                        # agregar puntos perpendiculaes 3D a la lista
                        lst_pntg = lst_pntg + [pntg_out1, pntg_out2]

    # escribir lista de puntos 3D a featureclass
    arcpy.CopyFeatures_management(lst_pntg, fc_out)
  
def getAngle(line, sr):
    '''Definir angulo de la linea'''
    first_pnt = line.firstPoint
    last_pnt = line.lastPoint
    pntg1 = arcpy.PointGeometry(first_pnt, sr)
    pntg2 = arcpy.PointGeometry(last_pnt, sr)
    tpl = pntg1.angleAndDistanceTo(pntg2, 'PLANAR')
    angle = tpl[0]
    return angle
  
if __name__ == '__main__':
    main()

 

There is a problem when you generate perpendicular points near the towers, when these points start to overlap and have different heights (see blue oval).

In that case you would need to delimit the generation of the points. This could be done by dividing the sections with the bisector of the angle:

This can also be established with some python (may post the tool in the future if there is any interest).

 

Once you have the correct 3D points, you can create a TIN from them and translate the TIN to raster.

 

3D Python python snippets

Outcomes