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

[RachelAlbritton]

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)

[MartinSchoner]

Hi,
same here. I calculate viewsheds for Wind power plants. I got several points with OFFSETA (180m : the hight of the rotor) and OFFSETB (2m : the hight of an Observer)
Also mys results are not very realistic.

I notice this: the results differ according to the position oft the point in the raster cell. LL corner vs. UR corner of the raster cell makes a huge difference.
My thoughts: a "lee effect" by rhe raster cell itself with the OFFSETA

According to your experience that the 2nd floor sees lesser than the 1st floor I presume that the viewshed tool didn't work propperly with the OFFSETA - with points as an Input!
I changed the input: a polyline between the points togehter with the same parameters (OFFSET ...).
The result was more the way I expected it.

Another way is to skip the OFFSETA Parameter. At your obs. point add the floor height to the bare earth elevation raster.

regards
-------------------------------------------------------
   Friedemann Wagner , Geograph, M.A.
   Regionalverband Südlicher Oberrhein,
   Reichsgrafenstra�?e 19, 79102 Freiburg
   Tel.: +49(0)761/ 70327-24
   -------------------------------------------------------

rralbritton;340310 wrote:
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.