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
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:
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) # 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) # 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, row, row x2, y2, z2 = row, row, row 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