Determining height of polygon to block view of point as seen from polyline feature?

Question

With ArcMap 10.5, I have the following data:

1. a polyline (trail) from which a 2m observer will be looking;
2. a polygon (mound of dirt, consider it flat on top) whose height can change;
3. the height of point behind the polygon, 10m; and,
4. a DEM of the surface.

I found Combining a DEM file with building vector files? for incorporating the height of the polygon into the DEM.

How high does the polygon need to be in order to block the view to the point for the observer who is on the polyline (i.e. at any point, or walking along the polyline)?

Assume the polygon is shaped such that it will hide the point from any angle of view on the polyline.

I have used the Visibility tool previously (see screenshot). It answers a different question: from the polyline what points on the DEM raster are visible? And of course, an offset can be set for the observer height (my 2m) and surface height (my 10m).

Is there a way to find the height of the polygon?

Answer 1

I have created a very ugly, but functioning script to do what you are after. There are likely more elegant ways to do this, and my script probably has lots of semantic issues that people may pick apart. I imagine we will need to work through some of it to get it to work for you, but give it a go. I hope that this is down the lines of what you were looking for.

import arcpy
from arcpy import env
import gc
from arcpy.sa import *
arcpy.CheckOutExtension("spatial")
arcpy.env.outputMFlag = "Disabled"

polyline = r'C:\TEMP\polyline.shp'
polygon = r'C:\TEMP\polygon.shp'
point = r'C:\TEMP\point.shp'
raster = r'\\silver\projects\SpatialData_Core\Alberta\DEM\DEM_1_20000\Athabasca\\Athabasca\ras\dem_30'
output_location = r'C:\TEMP'

observer_height = 2
point_height = 10

arcpy.env.workspace = output_location
arcpy.env.overwriteOutput = True
dsc = arcpy.Describe(raster)
arcpy.env.extent = dsc.Extent
cellsize = (arcpy.GetRasterProperties_management(raster, 'CELLSIZEX')).getOutput(0)

## Create an analysis window to reduce processing time of loop
arcpy.Near_analysis(point, polyline,'','','','PLANAR')
max_distance = 0
for f in arcpy.SearchCursor(point):
    if f.getValue('near_dist') > max_distance:
        max_distance = f.getValue('near_dist')
arcpy.Buffer_analysis(point,'extent_buffer',max_distance * 2,'','','ALL') # Times 2 to try and ensure all of the polyline is contained within the buffer
arcpy.gp.ExtractByMask_sa(raster, 'extent_buffer.shp', 'extent_dem')

## Create an feature class to create new rasters
arcpy.Union_analysis('extent_buffer.shp;' + polygon, 'union_polygon')
arcpy.AddField_management('union.shp', 'height', 'SHORT')

## Turn point to polygon for conversion to raster
arcpy.Buffer_analysis(point,'buffer',cellsize,'','','ALL')
arcpy.AddField_management('buffer.shp', 'height', 'SHORT')
arcpy.CalculateField_management('buffer.shp', 'height', point_height)
arcpy.Union_analysis('extent_buffer.shp;buffer.shp', 'union_point')
arcpy.PolygonToRaster_conversion('union_point.shp', 'height', 'pt_adjust','','',cellsize)

x = 0
y = 0
while x == 0:
    print 'Testing polygon height of {} meters'.format(str(y))
    arcpy.MakeFeatureLayer_management('union.shp', 'layer1.lyr')
    arcpy.SelectLayerByAttribute_management('layer1.lyr','NEW_SELECTION', '"FID_POLYGO" > -1')
    arcpy.CalculateField_management('layer1.lyr', 'height', y)
    arcpy.SelectLayerByAttribute_management('layer1.lyr','CLEAR_SELECTION')
    arcpy.PolygonToRaster_conversion('union.shp', 'height', 'dem_adjust','','',cellsize)
    new_dem = Raster('extent_dem') + Raster('pt_adjust') + Raster('dem_adjust')
    arcpy.gp.Visibility_sa(new_dem, polyline, 'visibility', '', 'FREQUENCY','ZERO', '','','','','', observer_height)
    arcpy.gp.ExtractValuesToPoints_sa(point, 'visibility', 'visi_check')
    for f in arcpy.SearchCursor('visi_check.shp'):
        if f.getValue('rastervalu') == 0:
            print 'The polygon height that will block view of the point is {} meters.'.format(str(y))
            x = 1
        else:
            x = 0
            break
    if y == 10: # bound the number of iterations the loop can test
        print 'Polygon is up to 10 meters high and still the point is visible. Either increase the number of iterations I can run, or reconsider the placement/length of your polyline.'
        sys.exit()
    y += 1