rralbritton

PLease Help - Viewshed - Why do my results seem incorrect???

Discussion created by rralbritton on Oct 29, 2013
Latest reply on Dec 11, 2013 by RVSO
I've written a script that calculates the oceanfront viewshed of a parcel for both the first and second floor. The script runs fine, however, the results seem incorrect - mainly that a second floor viewshed is significantly less then the first floor which seems really off. Am I doing something incorrectly in my script?

Note - there a two rasters being used in this script, a lidar raster that contains elevations of the buildings, and a lidar bare earth elevation raster that is being used to calculate the bare earth elevation height of each parcel.

Also - my script is to long to post the entire thing, so i've attached it to this post as well.

The general methodology:

1. Loop through the observation point layer and select each point one by one
2. Select the observation points building footprint that it intersects with and convert to a temporary raster.
3. Set the elevation of the parcel building footprint to zero and merge with a lidar raster (integer, grid).
4. Using a bare earth lidar elevation raster (interger, GRID format) use the summary statistics tool to calculate the average elevation of the building foot print.
5. Update the SPOT value in the observation point file with the average bare earth elevation of the building footprint.
6. Set the OFFSETA value to correct height (10 ft. for first floor - 20 ft. for second floor)
7. Calculate a temporary viewshed using a lidar elevation raster (integer, grid) that contains the merged temporary parcel raster.
8. Convert results to polygon and clip to 1 mile buffer.
9. Calculate viewshed in degrees.

import arcpy, os, shutil, datetime
from arcpy import env
from arcpy.sa import*

#set workspaces
arcpy.env.workspace = arcpy.GetParameterAsText(0) 
outputWorkspace = arcpy.GetParameterAsText(1) 
arcpy.env.overwriteOutput = True

#Variables
ObsPts = arcpy.GetParameterAsText(2)
footprint =  arcpy.GetParameterAsText(3) 
Elevation = arcpy.GetParameterAsText(4)
BareElevation = arcpy.GetParameterAsText(5)
Ocean = arcpy.GetParameterAsText(6)
FloorField = arcpy.GetParameterAsText(7) 
Year = arcpy.GetParameterAsText(8)

def Viewshed (LAYER, FC, HEIGHT, FLOOR, FOLDER):
    """Returns the viewshed, in degrees for each input observation point,
    for first (10ft) and second (20ft) floors."""
    #FOLDER = Intermediate files that can be deleted
    
    before = getTime()#Gets the time each viewshed analysis starts
    
    #Select observation points one by one
    arcpy.SetProgressorLabel("Starting viewshed analysis...")
    RangeCount = int(arcpy.GetCount_management(LAYER).getOutput(0))
    
    arcpy.CalculateField_management(LAYER,"OFFSETA",HEIGHT,"PYTHON_9.3")
    arcpy.AddMessage("Observation point heights have been set to "+str(HEIGHT)+" ft.")
    
    for points in range (0,RangeCount):

        try:
            
            arcpy.AddMessage("Calculating viewshed for "+str(RangeCount)+" points, for floor number "+str(FLOOR))
            
            FID = "FID=%s" % (points)

            #Get bare earth elevation of parcel
            arcpy.SelectLayerByAttribute_management(LAYER,"NEW_SELECTION",FID)
            arcpy.SelectLayerByLocation_management(footprintFL,"INTERSECT",LAYER)
            outExtractByMask = ExtractByMask(BareElevation,footprintFL)
            outExtractByMask.save(FOLDER+"\\tempfp_"+str(FID)[4:])
            ElevAvg = ElevAvgTables+"\\avgelev_"+str(FID)[4:]+".dbf"
            arcpy.Statistics_analysis(outExtractByMask,ElevAvg,[["VALUE","MEAN"]])#gets mean bare earth elevation value of parcel
            arcpy.AddField_management(ElevAvg,"Pt_FID","SHORT")
            arcpy.CalculateField_management(ElevAvg,"Pt_FID",FID[4:])
            arcpy.AddJoin_management(LAYER,"FID",ElevAvg,"Pt_FID","KEEP_COMMON")
            Field1 = os.path.basename(FC).split(".")[0]+".SPOT"
            Field2 = "!"+os.path.basename(ElevAvg).split(".")[0]+".MEAN_VALUE!"
            arcpy.CalculateField_management(LAYER,Field1,Field2,"PYTHON_9.3")#adds teh mean earth bare elevation value to SPOT field in point layer
            arcpy.RemoveJoin_management(LAYER)

            arcpy.env.workspace=arcpy.GetParameterAsText(0)
            
            #Set parcel elevation to 0
            arcpy.SelectLayerByAttribute_management(LAYER,"NEW_SELECTION",FID)
            arcpy.SelectLayerByLocation_management(footprintFL,"INTERSECT",LAYER)
            RastFootprint = FOLDER+"\\fp_"+str(FID)[4:].split(".")[0]
            arcpy.PolygonToRaster_conversion(footprintFL,"FID",RastFootprint,"MAXIMUM_AREA","",6)
            outIsNull = IsNull(RastFootprint) #Identify NoData Areas
            outIsNull.save(FOLDER+"\\"+str(FID)[4:]+"_null")
            outCon = Con(outIsNull,Elevation,0) #Use Con tool to change building footprint to elevation of 0 while leaving all other building footprints as is
            outCon.save(FOLDER+"\\"+str(FID)[4:]+"_con")
            
            #perform viewshed analysis
            arcpy.SetProgressorLabel("Changing elevation footprint to bare earth elevation for point "+str(FID))
            arcpy.SetProgressorLabel("Performing Viewshed Analysis for point "+str(FID))
            outViewshed = FOLDER+"/"+"view_"+str(FID)[4:].split(".")[0]
            arcpy.Viewshed_3d(outCon,LAYER,outViewshed,1)

            #convert viewshed to polygon
            arcpy.SetProgressorLabel("Converting viewshed"+str(FID)+" to polygon")
            OutPoly = FOLDER+"/"+os.path.basename(outViewshed).split(".")[0]+"_poly.shp"
            arcpy.RasterToPolygon_conversion(outViewshed,OutPoly)

            #buffer selected viewpoint
            arcpy.SetProgressorLabel("Buffering point "+str(FID))
            outBuffer = FOLDER+"/buffer_"+str(FID)[4:].split(".")[0]+".shp"
            arcpy.Buffer_analysis(LAYER,outBuffer,"1 mile")
            
            #Clip buffer to Ocean
            arcpy.SetProgressorLabel("Clipping point "+str(FID)+" buffer to ocean")
            BufferClip = FOLDER+"/buffer_clipped"+str(FID)[4:].split(".")[0]+".shp"
            arcpy.Clip_analysis(outBuffer, Ocean, BufferClip)
            
            #Intersect viewshed Poly to Clipped Observation Point Buffer
            arcpy.SetProgressorLabel("Intersecting polygon viewshed" +str(FID)+ "with clipped buffer")
            ViewshedClip = Final_Floor_Viewsheds+"/FinalView_Flr"+str(FLOOR)+"_"+str(FID)[4:].split(".")[0]+".shp"
            arcpy.Intersect_analysis([BufferClip,OutPoly],ViewshedClip)
            
            arcpy.AddField_management(ViewshedClip,"LENGTH","DOUBLE")
            arcpy.CalculateField_management(ViewshedClip,"LENGTH","!SHAPE.length@miles!","PYTHON_9.3")
            
            #Run Summary Statistics for each viewshed
            arcpy.SetProgressorLabel("Running summary statistics and calculating viewshed in degrees for "+str(FID))
            SumTab = SummaryTables+"/summary_"+str(FLOOR)+"_"+str(FID)[4:].split(".")[0]+".dbf"
            StatsField = [["LENGTH","MEAN"]]
            CaseField = ["GRIDCODE","ID","SPOT","OFFSETA",FloorField] #will need to make the id field user input
            arcpy.Statistics_analysis(ViewshedClip,SumTab,StatsField,CaseField)

            #Add Field and Caluclate Viewshed in degrees
            arcpy.AddField_management(SumTab,"FLR_RAN","SHORT")#will hold the floor number currently being analyzed - 1 or 2
            arcpy.AddField_management(SumTab,"VIEW_"+Year,"DOUBLE")
            arcpy.CalculateField_management(SumTab,"VIEW_"+Year,"((!MEAN_LENGT!/6.28)*180)","PYTHON_9.3")
            arcpy.CalculateField_management(SumTab,"FLR_RAN",FLOOR)
            
            uc = arcpy.UpdateCursor(SumTab)
            View = "VIEW_"+Year
            for line in uc:
                if line.GRIDCODE ==0:
                    line.setValue(View,0)
                    uc.updateRow(line)
                del line
            del uc

            count = int(arcpy.GetCount_management(SumTab).getOutput(0))

            if count == 2:
                uc=arcpy.UpdateCursor (SumTab)
                for line in uc:
                    if line.GRIDCODE ==0:
                        uc.deleteRow(line)
                        break
                    del line
                del uc

            arcpy.SelectLayerByAttribute_management(LAYER,"CLEAR_SELECTION")
            arcpy.SelectLayerByAttribute_management(footprintFL,"CLEAR_SELECTION")
            arcpy.AddMessage("The viewshed for point "+str(FID)[4:]+" of "+str(RangeCount)+" on floor number "+str(FLOOR)+" is complete")

            after = getTime()#get time viewshed completed for individual point
            diffTime(before,after)

        except:

            arcpy.AddMessage("***An error occured while processing observation point "+str(FID)+"\nPlease refer to the error log located in "+outputWorkspace+" for details.***\n")
            arcpy.SelectLayerByAttribute_management(LAYER,"CLEAR_SELECTION")
        

TotalBefore = getTime()
    
#Create Featue Layers
arcpy.SetProgressorLabel("Creating feature layers")
PointsFL = outName(ObsPts,"_Lyr")
footprintFL = outName(footprint,"_Lyr")
arcpy.MakeFeatureLayer_management(ObsPts, PointsFL)
arcpy.MakeFeatureLayer_management(footprint,footprintFL)
   
Viewshed (PointsFL, ObsPts, 10, 1, IntermediateFiles)

Outcomes