Site viewshed with DSM: removing "false" visibility readings from trees and buildings

Question

I've been working on modeling a proposed building and running a visibility analysis to determine from what areas the building could potentially be seen using ArcGIS Desktop with 3D Analyst. Using a DSM raster derived from LiDAR (includes all points - ground, trees, buildings, etc.) as the base surface information so that above-ground sight obstructions are accounted for.

My process to find the areas that the new building could be seen from has been to set the proposed building roof corners and peaks as 3D points and use them as "observers" - if the roof corner can "see" an area, then the building can be seen from that area. I used the Visibility tool in the 3D Analyst toolbox and Model Builder to run through the observer points at the building roof (16 all together) to generate visibility rasters and then mosaic them together.

The issue is that using the DSM is a double-edged sword - using it to account for above-ground obstructions has also created speckled areas of "visibility" on treetops and building tops where the new building could be seen from if the observer's eyes were at that height. But obviously, an observer would not be standing on a rooftop or treetop, so these areas are erroneous. See image below (hillshade with red areas being visible from the proposed building site)

How can I remove these areas, besides manual cleanup? The LiDAR is not classified for vegetation and buildings, only ground and water, so using those LiDAR points as a mask to delete the erroneous areas would require some classification work.

Answer 1

This was something I'd been thinking on for viewsheds, didn't realize Line of Sight was going to need z values. I included a global timer to find out if there is a particular step that slows it down. Its written to reduce by an offset factor, but I included comments on how to alter it to use variable heights from another raster instead, like a DEM.

# Version: ArcGIS 10.2
import sys, os
import arcpy
from arcpy import env
from arcpy import da
from arcpy.sa import *

arcpy.CheckOutExtension("3D")
arcpy.CheckOutExtension("Spatial")
arcpy.env.parallelProcessingFactor = "100%"

#Global Timer
def exec_time(start, message):
    end = time.clock()
    comp_time = end - start
    arcpy.AddMessage("Run time for " + message + ": " + str(comp_time))
    start = time.clock()
    return start
###Please set these params###
#set your workspace
path=
env.workspace = path
#input Elev Raster
elevRaster =
#input Raster (viewshed result)
visibleRaster =
#input observer points
observerPnts =
#intermediate observer points with z values
observerPntsZ=
#intermediate visible points, and with z values
visPnts =
visPntsZ =
#output sight lines
visLn =
#ouput line of sight lines
resultsLn =
#Use adjFactor as correction height, adjRaster could be used instead for variable  adjustment factor from a DEM
adjFactor =
######
start=exec_time(start, "startup")
#convert the raster into set of points for all visible locations
OutRasVis = SetNull(visibleRaster, visibleRaster, "VALUE = 0") 
OutRasVis.save(path + "\\tempRasVis")
start=exec_time(start, "removing nulls")
arcpy.RasterToPoint_conversion(OutRasVis, visPnts)
start=exec_time(start, "created visible points")
#set heights for observers based on elevation raster (wasn't sure where this info was)
arcpy.MakeFeatureLayer_management(observerPnts, observerPntsLyr)
ExtractValuesToPoints(observerPntsLyr, elevRaster, observerPntsZ, "NONE", "VALUE_ONLY") 
#set heights for visible, use adjRaster in place of elevRaster if location based adjustment
arcpy.MakeFeatureLayer_management(visPnts, visPntsLyr)
ExtractValuesToPoints(visPntsLyr, elevRaster, visPntsZ, "NONE", "VALUE_ONLY")
start=exec_time(start, "elevations to points")

###skip this if using location based adjustment###
arcpy.AddField_management(visPntsZ, "Spot", "Double")
with arcpy.da.UpdateCursor(visPntsZ, ["RASTERVALU", "Spot"]) as cursor:
    for row in cursor:
        if (row[1] == 0):
            row[1]= row[0] - adjFactor
            cursor.updateRow(row)
start=exec_time(start, "adjustment factor")
###skip these is using location based adjustment###

#create observer -> visible sight lines, change "Spot" to "RASTERVALU" if using location based adjustment
arcpy.ConstructSightLines_3d (observerPntsZ, visPntsZ, visLn, "RASTERVALU", "Spot","","","")
arcpy.LineOfSight_3d(elevRaster, outLn, resultsLn, "","","","","")
#delete undesired intermediates (visPnts, visPntsZ, observerPntsZ, visLN)
arcpy.Delete_management(visPnts)

I've not debugged it, so if you find anything let me know and I'll update this code for anyone in the future. One concern I have about this approach is that since the viewer location is lowered but not the surrounding cells any interior cells will no longer be visible(i.e. only the outside perimeter of a house can see out) vastly reducing viewers.

Answer 2

False positives in Viewsheds is something i have tried to highlight a few times including in a GISRUK 2012 poster (available on my Linkedin page - David Roberts-Lock), to resolve it i always obtain a LiDAR DTM model as well and then (using Raster Calculator on the two esri grids do DSM - DTM) this gives me a raster i class as Obstructions (or surface features if you will). once obtained (and depending on the quality of the LiDAR survey) you reclassify to find just areas where the difference is greater than say 0.5 meters. making all the values above 0.5 noData and those below a 0 (zero). adding this raster to the Viewshed result will then remove the features but ensure they are still hiding visibility from the terrain behind.

This all needs Spatial Analyst obviously but means coding isnt required!! (although you can make into an Esri Tool using Modelbuilder if you really want)