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 #Check out the ArcGIS 3D Analyst extension license arcpy.CheckOutExtension("3D") arcpy.CheckOutExtension("Spatial") arcpy.SetProgressor("default","Conducting Viewshed Analysis") #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) #Set analysis extent to elevation raster arcpy.env.extent = Elevation arcpy.env.cellSize = Elevation #Open error log file infile = open(outputWorkspace+"\\Error_Log_"+Year+".txt","w") #Count number of parcels being processed arcpy.AddMessage("\nCalculating viewshed for "+str(RangeCount)+" parcels") sc = arcpy.SearchCursor(PointsFL) for row in sc: try: #Get Parcel ID value value = row.ID count = row.FID+1 FlrCnt = row.getValue(FloorField) #Get bare earth elevation of parcel arcpy.SetProgressorLabel("Changing elevation footprint to bare earth elevation for point "+str(value)) SQL = "Id =" +str(value)+"" arcpy.SelectLayerByAttribute_management(PointsFL,"NEW_SELECTION",SQL) arcpy.SelectLayerByLocation_management(footprintFL,"INTERSECT",PointsFL) arcpy.env.workspace = IntermediateFiles #need to change workspace so that the .save files get saved correctly outExtractByMask = ExtractByMask(BareElevation,footprintFL) outExtractByMask.save(IntermediateFiles+"\\ebm_"+str(value)) ElevAvg = ElevAvgTables+"\\avgelev_"+str(value)+".dbf" arcpy.Statistics_analysis(outExtractByMask,ElevAvg,[["VALUE","MEAN"]]) arcpy.AddField_management(ElevAvg,"Pt_ID","SHORT") arcpy.CalculateField_management(ElevAvg,"Pt_ID",value) arcpy.AddJoin_management(PointsFL,"Id",ElevAvg,"Pt_ID","KEEP_COMMON") Field1 = os.path.basename(ObsPts).split(".")[0]+".SPOT" Field2 = "!"+os.path.basename(ElevAvg).split(".")[0]+".MEAN_VALUE!" arcpy.CalculateField_management(PointsFL,Field1,Field2,"PYTHON_9.3") arcpy.RemoveJoin_management(PointsFL) #Set parcel elevation to 0 this will be replaced by SPOT value caclualted above RastFootprint = IntermediateFiles+"\\fp_"+str(value).split(".")[0] arcpy.PolygonToRaster_conversion(footprintFL,"FID",RastFootprint,"MAXIMUM_AREA","",6) outIsNull = IsNull(RastFootprint) #Identify NoData Areas outIsNull.save(IntermediateFiles+"\\null1_"+str(value)) 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(IntermediateFiles+"\\con1_"+str(value)) #Final raster to be used in viewshed analysis #buffer selected viewpoint arcpy.SetProgressorLabel("Buffering point "+str(value)) outBuffer = IntermediateFiles+"\\buffer_"+str(value)+".shp" arcpy.Buffer_analysis(PointsFL,outBuffer,"1 mile") #Convert buffer polygon to line BufferLine = IntermediateFiles+"\\BufferLine_"+str(value)+".shp" arcpy.FeatureToLine_management(outBuffer,BufferLine) #Clip buffer to Ocean arcpy.SetProgressorLabel("Clipping point "+str(value)+" buffer to ocean") BufferClip = IntermediateFiles+"\\buffer_clipped_"+str(value).split(".")[0]+".shp" arcpy.Clip_analysis(outBuffer, Ocean, BufferClip) if FlrCnt ==1: #parcel floor count =1 arcpy.AddMessage("\nParcel "+str(value)+" has 1 story to process. Calculating viewshed now...") print "\nParcel ",str(value)," has 1 story to process. Calculating viewshed now..." DegViewshed(1,10) #Calculate the viewshed with an observer height of 10 feet then move to point arcpy.AddMessage("First floor viewshed for parcel "+str(value)+" has been completed...") print "First floor viewshed for parcel ",str(value)," has been completed..." arcpy.AddMessage(str(count)+" of "+str(RangeCount)+"parcles has been completed.\n") print str(count)," of "+str(RangeCount)," parcels has been processed.\n" else: #if parcel has 1.5 floors or greater do this arcpy.AddMessage("\nParcel "+str(value)+" has 2 stories to process. Calculating viewsheds now...") print "\nParcel ",str(value)," has 2 stories to process. Calculating viewsheds now..." DegViewshed(1,10) #Calculate first floor view, then arcpy.AddMessage("First floor viewshed for parcel "+str(value)+" has been completed...") print "First floor viewshed for parcel ",str(value)," has been completed..." DegViewshed(2,20) #Calculate second floor view arcpy.AddMessage("Second floor viewshed for parcel "+str(value)+" has been completed...") print "Second floor viewshed for parcel ",str(value)," has been completed..." arcpy.AddMessage("Viewsheds for "+str(count)+" of "+str(RangeCount)+" parcels have been processed.\n") print "Viewsheds for",str(count)," of ",str(RangeCount),"parcels have been processed.\n" except: arcpy.AddMessage("***An error occured processing parcel "+str(value)+". Refer to error log for details.") print "***An error occured processing parcel "+str(value)+". Refer to error log for details." infile.write("An error occured processing parcel "+str(value)+".\n") infile.write(arcpy.GetMessages()+"\n") arcpy.SelectLayerByAttribute_management(PointsFL, "CLEAR_SELECTION") arcpy.SelectLayerByLocation_management(footprintFL,"CLEAR_SELECTION") del row del sc del row del sc
Solved! Go to Solution.
except: pass
except: arcpy.AddMessage("***An error occured processing parcel "+str(value)+". Refer to error log for details.") print "***An error occured processing parcel "+str(value)+". Refer to error log for details." infile.write("An error occured processing parcel "+str(value)+".\n") infile.write(arcpy.GetMessages()+"\n") arcpy.SelectLayerByAttribute_management(PointsFL, "CLEAR_SELECTION") arcpy.SelectLayerByLocation_management(footprintFL,"CLEAR_SELECTION") break del row del sc continue
except: pass