As mentioned in my previous post Working with 3D and M-aware geometries in Arcpy, you can create 3D features with Arcpy. We needed those for a project that involved marking trees for logging that interfere with the design of overhead power lines.  One of the conditions is that the tree should not be closer than 6 meters to the cables. In the 3D Analyst toolbox there is a tool “Near 3D” that calculates the 3D distance between 3D features. Exactly what we need to solve this type of problem. The challenge however, is to create 3D features to be able to use this tool.

The forest inventory yields a point featureclass with all kinds of data on the tree. Additional research provides potential growth on crown and height. Based on these specifications we can create 3D features of the trees that can be used in combinations with the 3D overhead power lines to determine the 3D distance between the potential size of the trees and the lines.

In case we use a very simple model of a tree, it consists of the trunk and the crown. We will need the following data to start with.

• Diameter of trunk
• Height of trunk
• Diameter  crown
• Height of crown
• Elevation

So, let’s start with a little python code.  The idea is to create rings (polylines) for the trunk and store these as a multipart in the output featureclass. For the crown we will also create rings and use a sine function to simulate the size of the diameter at each elevation step. The result will be something like this:

Observe the code below:

import arcpy

import math

def main():
    import os

    # output
    fc_out = r'C:\GeoNet\_WorkingWith3D\data.gdb\Tree3D_v01'
    arcpy.env.overwriteOutput = True

    # the tree point
    sr = arcpy.SpatialReference(3857)
    x = -13046027.327
    y = 4036381.962
    elevation = 365

    # dimensions of tree
    tree_height = 30
    crown_diameter = 20
    trunk_diameter = 0.7
    trunk_height = 18
    crown_height = tree_height - trunk_height

    # configuration of number of points and lines (detail of tree)
    crown_lines = 20  # number of lines for crown
    trunk_lines = 25 # number of lines for trunk
    pnts_on_circle = 18  # a point every 20 degrees on circle

    # create output featureclass
    ws, fc_name = os.path.split(fc_out)
    arcpy.CreateFeatureclass_management(ws, fc_name, "POLYLINE", None, None, "ENABLED", sr)

    # Addional field for tree part to store Trunk or Crown
    fld_tree_part = "TreePart"
    arcpy.AddField_management(fc_out, fld_tree_part, "TEXT", None, None, 10)

    # insert cursor to write feature
    flds_out = ('SHAPE@', fld_tree_part)
    with arcpy.da.InsertCursor(fc_out, flds_out) as curs_out:

        # list to store multipart feature
        lst_feats = []
        pnt = arcpy.Point(x, y, elevation)
        pntg = arcpy.PointGeometry(pnt, sr, True)

        # define z values
        z0 = elevation
        z1 = z0 + trunk_height
        z2 = z1 + crown_height

        # create multipart feature for trunk
        polyline = createTree3Dtrunk(pntg, trunk_lines, pnts_on_circle, z0, z1, trunk_diameter, sr)
        # write to output
        row_out = (polyline, "Trunk", )
        curs_out.insertRow(row_out)

        # create multipart feature for crown
        polyline = createTree3Dcrown(pntg, crown_lines, pnts_on_circle, z1, z2, crown_diameter, sr)
        # write to output
        row_out = (polyline, "Crown", )
        curs_out.insertRow(row_out)


def createTree3Dtrunk(pntg, trunk_lines, pnts_on_circle, z0, z1, trunk_diameter, sr):
    ''''create line multipart feature
        each line represents a circle of the trunk'''
    lst_lines = []
    for i in range(trunk_lines):
        # determine z value of trunk line
        zi = float(i) / float(trunk_lines - 1) * (z1  - z0) + z0
        # create line by buffering
        circ = pntg.buffer(trunk_diameter / 2.0)
        # get outline (boundary) of buffer polygon
        circ_outline = circ.boundary()

        lst_pnts = []
        # create points on circle
        for a in range(pnts_on_circle):
            # determine fraction of line for pnt
            f = float(a) / float(pnts_on_circle)
            pnt_on_line = circ_outline.positionAlongLine(f, True)
            # insert the z value for this trunk line
            pnt_on_line_z = arcpy.Point(pnt_on_line.firstPoint.X, pnt_on_line.firstPoint.Y, zi)
            # append point to list
            lst_pnts.append(pnt_on_line_z)

        # add first point to close line
        lst_pnts.append(lst_pnts[0])
        # add list of points to list of lines (for multipart)
        lst_lines.append(lst_pnts)

    # create multipart feature
    polyline = arcpy.Polyline(arcpy.Array(lst_lines), sr, True, False)
    return polyline


def createTree3Dcrown(pntg, crown_lines, pnts_on_circle, z1, z2, crown_diameter, sr):
    ''''create line multipart feature
        each line represents a circle of the crown'''
    lst_lines = []
    for i in range(1, crown_lines):
        print i
        # determine z value of crown line
        zi = float(i) / float(crown_lines) * (z2  - z1) + z1
        print zi
        # determine angle expressed in pi to determine diameter of crown at zi
        angle_pi = float(i) / (crown_lines) * math.pi
        print angle_pi
        diam_i = math.sin(angle_pi) * float(crown_diameter)
        print angle_pi
        # create line by buffering
        circ = pntg.buffer(diam_i / 2.0)
        # get outline (boundary) of buffer polygon
        circ_outline = circ.boundary()

        lst_pnts = []
        # create points on circle
        for a in range(pnts_on_circle):
            # determine fraction of line for pnt
            f = float(a) / float(pnts_on_circle)
            print " - f", f
            pnt_on_line = circ_outline.positionAlongLine(f, True)
            # insert the z value for this trunk line
            pnt_on_line_z = arcpy.Point(pnt_on_line.firstPoint.X, pnt_on_line.firstPoint.Y, zi)
            # append point to list
            lst_pnts.append(pnt_on_line_z)

        # add first point to close line
        lst_pnts.append(lst_pnts[0])
        # add list of points to list of lines (for multipart)
        lst_lines.append(lst_pnts)

    # create multipart feature
    polyline = arcpy.Polyline(arcpy.Array(lst_lines), sr, True, False)
    return polyline

if __name__ == '__main__':
    main()

There are two functions that do the work: createTree3Dtrunk  and createTree3Dcrown. The first one creates the multipart trunk polyline and the second one the multipart crown polyline.

The first part of the script defines the dimensions of the tree and creates the empty output featureclass (Z enabled) and adds an additional field that will the text "Trunk" or "Crown".  After defining the insert cursor, the point geometry is created and the trunk and crown heights are translated to elevation.

To create the trunk lines the createTree3Dtrunk  function is called. This function iterates over the number of trunk lines, and determines the elevation for that line. A buffer is calculated to create the circle. You may wonder why this line is not used directly as part in the multipart polyline, but creating a buffer based on a 3D point geometry will return a 2D polygon. Therefore, points are located on the boundary of the buffer, and the point is upgraded to 3D and a 3D polyline is created with these points. Each list of points for a single line is appended  to a list of lists, and this is the list which is used to create the multipart 3D polyline.

For the crown a similar method is used, however, the size of the buffer is defined by the sine of the fraction of elevation (between where the crown starts and end). That part is translated into 0 - 1π creating the proper effect of the change in size of the crown.

This geometry (or better yet, a geometry containing both the trunk and the crown) can be used in the Near 3D (3D Analyst Toolbox) tool to determine the smallest distance between any part of the tree and the nearest 3D overhead power line.

The 3D distance can be used to visualize the 3D distance between the tree and the overhead power line:

The optional location will yield additional fields to indicate which location of the tree is nearest to the closest location on the overhead power line. These fields can be used to generate lines between those two points in case you want to explore the distances in more detail. Below a snippet of code generate these 3D lines:

def main():
    import arcpy

    fc =r'D:\Test\data.gdb\TreePolygon3D_v01'
    fc_out = r'D:\Test\data.gdb\near3Dlines_v01'
    sr = arcpy.Describe(fc).spatialReference

    i = 0
    lst_feats = []
    flds = ('NEAR_FROMX', 'NEAR_FROMY', 'NEAR_FROMZ',
            'NEAR_X', 'NEAR_Y', 'NEAR_Z')
    with arcpy.da.SearchCursor(fc, flds) as curs:
        for row in curs:
            i += 1
            if i % 100 == 0:
                print "Processing:", i
            x1, y1, z1 = row[0], row[1], row[2]
            x2, y2, z2 = row[3], row[4], row[5]

            pnt1 = arcpy.Point(x1, y1, z1)
            pnt2 = arcpy.Point(x2, y2, z2)
            polyline = arcpy.Polyline(arcpy.Array([pnt1, pnt2]), sr, True, False)
            lst_feats.append(polyline)

    arcpy.CopyFeatures_management(lst_feats, fc_out)


if __name__ == '__main__':
    main()

Optionally calculate the 3D length of the lines and use that to visualize the lines in different length classes:

Next up: create a tree based on polygons

Python python snippets 3D