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) vis

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.

  • Thoughts on taking your set of locations that are visible from the building and running line of site from those points to the building points with an adjusted elevation for the visible point? Depending on the scope of your area it could get intensive fast, are you open to a python solution? – jbosq Feb 13 '15 at 19:45
  • The visible areas cover about a 1 mile radius, so it's not too large, but not small either. Could perhaps create gridded points within the visible areas, offset them 5' from the DEM for eye level, and batch run line of sight, but like you said, that could be some major processing time. I'm open to python – mpianka Feb 13 '15 at 20:21

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.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
env.workspace = path
#input Elev Raster
elevRaster =
#input Raster (viewshed result)
visibleRaster =
#input observer points
observerPnts =
#intermediate observer points with z values
#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
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)

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.


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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.