Python code for buffer creation and intersection

Question

The script below uses an ArcPy package which helps to create buffer of the impact around a location where a blast occurs(user entered), intersects with building footprint which is the point of interest to figure out if the blast impact will be there on the point of interest. This I am intending to use for blast analysis and the impact it does on the point of interest.

So far the code runs fine without any issues, however, I was curious if someone could sense check the code and review it for a better way to write the same and any modifications that could help. It's a simple though a long code so I am not expecting anyone to rewrite the entire code, but any alternative style or corrections in error handling to make it more efficient will be helpful.

import arcpy
import sys,os
import datetime
from datetime import datetime
arcpy.env.overwriteOutput = True
##Tool input##
Pfc = arcpy.GetParameterAsText(0)##Blast locations
idField = arcpy.GetParameterAsText(1) ##Blast ID field - can be object id or user defined/named field
pointID = arcpy.GetParameterAsText(2) ##ID of blast location to be used (Amazon = Front Entrance or Loading Bay)
bFootprints = arcpy.GetParameterAsText(3)##Building footprint
totalEx = arcpy.GetParameterAsText(4) ##Total building exposure
buildType = arcpy.GetParameterAsText(5) ##Building construction type - this influences loss ratios
outGDB = arcpy.GetParameterAsText(6)##Output Location (gdb)
buffName = arcpy.GetParameterAsText(7) ##Buffer name - this will persist for intersect name(s) too
BlastZones = arcpy.GetParameterAsText(8) ##Blast zones to use - DHS or Lloyds.
BlastSize = arcpy.GetParameterAsText(9) ##Blast size - drop down list
outCorName = arcpy.GetParameterAsText(10) ##output coordinate system - drop down list
lossTable = arcpy.GetParameterAsText(11) ##Output loss table (to GDB)
outShape = arcpy.GetParameterAsText(12) ##Output to shapefile - T or F
outCSV = arcpy.GetParameterAsText(13) ##Output loss table as csv - T or F
##NAD 83
NAD83 = arcpy.SpatialReference(6350)
##WGS 84 Web Mercator
WebMerc = arcpy.SpatialReference(3857)
##BNG
BNG = arcpy.SpatialReference(27700)
##Check if being run in ArcMap - if yes, change some variables
if 'ArcMap.exe' in os.path.basename(sys.executable):
 ##Alter Pfc
 desc = arcpy.Describe(Pfc)
 Pfc = str(desc.path) + "\\" + desc.name
 ##Alter fc
 desc = arcpy.Describe(Pfc)
 fc = str(desc.path) + "\\" + desc.name
####Processing starts here...####
##Set output coordinates
if outCorName == "NAD83":
 outCor = NAD83
elif outCorName == "WebMerc":
 outCor = WebMerc
elif outCorName == "BNG":
 outCor = BNG
##Set csv and shapefile output location
outFolder = os.path.join(outGDB.rsplit("\\",1)[0],str(buffName+"_output"))
if os.path.exists(outFolder):
 pass
else:
 arcpy.CreateFolder_management(outGDB.rsplit("\\",1)[0],str(buffName+"_output"))
##UID field - used in update cursor
BF_name = bFootprints.rsplit("\\",1)[1]
####Set blast size(s) and associated loss ratios
if BlastZones == "Lloyds":
 if BlastSize == "Backpack":
 blastSize = "36;65;104"
 LR1 = 0.8
 LR2 = 0.4
 LR3 = 0.05
 elif BlastSize == "truck1Ton":
 blastSize = "121;213;338"
 LR1 = 0.8
 LR2 = 0.4
 LR3 = 0.05
 elif BlastSize == "truck2Ton":
 blastSize = "656;1312;1640"
 LR1 = 1
 LR2 = 0.25
 LR3 = 0.1
if BlastZones == "DHS":
 if BlastSize == "Backpack":
 blastSize = "26;80;125"
 LR1 = 1
 LR2 = 0.25
 LR3 = 0.1
 elif BlastSize == "truck1Ton":
 blastSize = "85;260;375"
 LR1 = 1
 LR2 = 0.25
 LR3 = 0.1
 elif BlastSize == "truck2Ton":
 blastSize = "100;335;500"
 LR1 = 1
 LR2 = 0.25
 LR3 = 0.1
 elif BlastSize == "truck5Ton":
 blastSize = "140;450;640"
 LR1 = 1
 LR2 = 0.25
 LR3 = 0.1
 elif BlastSize == "truck6Ton":
 blastSize = "175;590;980"
 LR1 = 1
 LR2 = 0.25
 LR3 = 0.1
##Set for area calc.
geom = "AREA_GEODESIC"
##Add new fields to building footprints - to store intersected blast radius
tmp = blastSize.rsplit(";")
Field1 = "BuildingArea"
Field2 = "AREA_GEO"
fList = [Field1,Field2]
##sList = [Field1,Field2,Field3,Field4]
TField1 = "DamageLevel"
TField2 = "ImpactZoneByFeet"
TField3 = "PropertyFootprint_Area"
TField4 = "PropertyImpacted_Area"
TField5 = "TotalExposure"
TField6 = "ImpactRatio"
TField7 = "ImpactExposure"
TField8 = "LossRatio"
TField9 = "Loss"
TablefList = [TField1,TField2,TField3,TField4,TField5,TField6,TField7,TField8,TField9]
tExp = float(totalEx)
IZ1 = tmp[0]
IZ2 = tmp[1]
IZ3 = tmp[2]
tableTemplate = r"O:\LossTable"
##Building intersection function##
def buildIntersect(outInt,fields,fc):
 with arcpy.da.SearchCursor(outInt,fields) as sCursor:
 for row in sCursor:
 bID = str(row[0])
 BR = int(row[1])
 intArea = row[2]
 where = "OBJECTID = " + bID
 with arcpy.da.UpdateCursor(fc,fList,where) as uCursor:
 for row in uCursor:
 if BR == int(tmp[0]):
 row[0] = intArea
 elif BR == int(tmp[1]):
 row[1] = intArea
 elif BR == int(tmp[2]):
 row[2] = intArea
 uCursor.updateRow(row)
 del sCursor, uCursor
##Loss Table function##
def LossTable(fc, TablefList, tableName):
 bArea = []
 intArea = []
 with arcpy.da.SearchCursor(fc,fList) as sCursor:
 for row in sCursor:
 bArea.append(row[0])
 intArea.append(row[1])
 del sCursor
 table = os.path.join(outGDB,tableName)
 arcpy.CopyRows_management(tableTemplate,table)
 with arcpy.da.UpdateCursor(table,TablefList) as uCursor:
 for row in uCursor:
 DL = row[0]
 if DL == "Heavy":
 row[1] = IZ1
 row[2] = bArea[0]
 row[3] = intArea[0]
 row[4] = totalEx
 row[5] = row[3]/row[2]
 row[6] = int(row[4])*row[5]
 row[7] = LR1
 row[8] = row[6]*row[7]
 elif DL == "Moderate":
 row[1] = IZ2
 row[2] = bArea[0]
 row[3] = intArea[1]
 row[4] = totalEx
 row[5] = row[3]/row[2]
 row[6] = int(row[4])*row[5]
 row[7] = LR2
 row[8] = row[6]*row[7]
 elif DL == "Light":
 row[1] = IZ3
 row[2] = bArea[0]
 row[3] = intArea[2]
 row[4] = totalEx
 row[5] = row[3]/row[2]
 row[6] = int(row[4])*row[5]
 row[7] = LR3
 row[8] = row[6]*row[7]
 uCursor.updateRow(row)
 del uCursor
 OutTableName = str(tableName) +".csv"
 if str(outCSV) == "true":
 arcpy.TableToTable_conversion(table,outFolder,OutTableName)
 arcpy.AddMessage("Loss calculation table created and saved in: " + outFolder)
##Output shapefile function##
def OutShapes(fcList,outFolder):
 for fc in fcList:
 outShp = str(fc.rsplit("\\",1)[1])+".shp"
 print(outShp)
 arcpy.FeatureClassToFeatureClass_conversion(fc,outFolder,outShp)
try:
 ##Blast Pressures
 Pres_1 = 10
 Pres_2 = 2
 Pres_3 = 1
 ####Processing starts here#### 
 fcList = []
 if type(pointID) == str:
 pointID = "'"+pointID+"'"
 print(pointID)
 field = str(idField)
 where = field +" = " + str(pointID) 
 memoryFeature = "in_memory" + "\\" + "myMemoryFeature"
 arcpy.Select_analysis(Pfc, memoryFeature, where)
 baseName = buffName
 buffName = os.path.join(outGDB,str(baseName + "_Buffer"))
 fcList.append(buffName)
 ##Add geometry in selected CRS to building footprints
 arcpy.AddGeometryAttributes_management(bFootprints,geom,"","",outCor)
 arcpy.AddMessage("Geometry added")
 ##Add field if needed
 if len(arcpy.ListFields(bFootprints,"BuildingArea"))>0:
 print("Building area field exists")
 else:
 arcpy.AddField_management(bFootprints, "BuildingArea", "DOUBLE")
 arcpy.CalculateField_management(bFootprints, "BuildingArea", "!AREA_GEO!",expression_type="PYTHON_9.3")
 fcList.append(bFootprints)
 print("Buffering...")
 arcpy.MultipleRingBuffer_analysis(memoryFeature, buffName, Distances=blastSize, Buffer_Unit="Feet", Field_Name="BlastRadius", Dissolve_Option="ALL", Outside_Polygons_Only="FULL")
 arcpy.AddMessage("Buffering done")
 print("Buffering done")
 tableName = baseName + "_LossTable"
 arcpy.Delete_management(memoryFeature)
 ##Intersect buffer with building footprint
 Buff = os.path.join(outGDB,str(baseName + "_Buffer"))
 outInt = os.path.join(outGDB,str(baseName+"_BuildingIntersect"))
 fcList.append(outInt)
 print("Intersecting...")
 arcpy.Intersect_analysis(in_features = [Buff, bFootprints], out_feature_class=outInt, join_attributes="ALL", output_type = "INPUT")
 arcpy.AddGeometryAttributes_management(outInt,geom,"","",outCor)
 ##Update bFootprints with intersected values where building ID and blast radius match
 arcpy.AddMessage("Intersecting complete")
 print("Intersecting complete")
 if str(lossTable) == "true":
 tableName = str(baseName) + "_LossTable"
 LossTable(outInt,TablefList,tableName)
 if str(outShape) == "true": 
 OutShapes(fcList,outFolder)
 arcpy.AddMessage("Feature classes copied to shapefile in the folder: " + outFolder)
except Exception, e:
 import traceback, sys, os
 tb = sys.exc_info()[2]
 print("Line %i" % tb.tb_lineno)
 arcpy.AddMessage("Line %i" % tb.tb_lineno)
 arcpy.AddMessage(str(e))
 arcpy.AddMessage(arcpy.GetMessages())
 print(str(e))
finally: 
 del arcpy

• Simon Forsberg: If you want to have this question deleted, the process for asking for deletion of an answered question is described in meta.stackexchange.com/a/5222/213556

Answer 1

Congratulations on the development. I have been able to verify some lines of code that can be better managed with the module "os" in order to take advantage of it more since you have it imported. I also suggest using the MakeFeatureLayer tool instead of Select_analysis, since the latter has a much longer processing time, in addition, the purpose in that section is to obtain a temporary layer in the processing. I also suggest assigning variables to the operations, to facilitate the use during the process. I suspect that in the section where the arcpy.da cursors are used, you can also optimize but for this it would be necessary to perform the tests with the ScripTool and real information. Excuse my English, I'm from Perú, greetings.

outFolder = os.path.abspath(os.path.join(outGDB, '../%s_output' % buffName))
# outFolder = os.path.join(outGDB.rsplit("\\",1)[0],str(buffName+"_output"))
#------------------------------------------------------------------------------
if not os.path.exists(outFolder):
 arcpy.CreateFolder_management(os.path.dirname(outFolder), os.path.basename(outFolder))
# if os.path.exists(outFolder):
# pass
# else:
# arcpy.CreateFolder_management(outGDB.rsplit("\\",1)[0],str(buffName+"_output"))
#------------------------------------------------------------------------------
##UID field - used in update cursor
BF_name = os.path.basename(bFootprints)
# BF_name = bFootprints.rsplit("\\",1)[1]
#------------------------------------------------------------------------------
print("Buffering...")
arcpy.MultipleRingBuffer_analysis(mfl, buffName, Distances=blastSize, 
Buffer_Unit="Feet", Field_Name="BlastRadius", Dissolve_Option="ALL", 
Outside_Polygons_Only="FULL")
# arcpy.MultipleRingBuffer_analysis(memoryFeature, buffName, Distances=blastSize, 
# Buffer_Unit="Feet", Field_Name="BlastRadius", Dissolve_Option="ALL", 
# Outside_Polygons_Only="FULL")
#------------------------------------------------------------------------------
arcpy.Delete_management(mfl)
# arcpy.Delete_management(memoryFeature)

Answer 2

Don't confuse module and module symbol imports

You wrote this:

import datetime
from datetime import datetime

which is both confusing and unnecessary. It doesn't seem that you use datetime anywhere, so just delete these lines. Even if you did need to use datetime, choose one or the other - either just import the module and use it for qualified names, or import symbols from it - not both.

Make some functions

You don't have any functions. Everything is hanging out in the same global scope. This is cluttered, difficult to read and difficult to maintain and debug. Write some functions!

Make the computer do the repetition

This:

Pfc = arcpy.GetParameterAsText(0)##Blast locations
idField = arcpy.GetParameterAsText(1) ##Blast ID field - can be object id or user defined/named field
pointID = arcpy.GetParameterAsText(2) ##ID of blast location to be used (Amazon = Front Entrance or Loading Bay)
bFootprints = arcpy.GetParameterAsText(3)##Building footprint
totalEx = arcpy.GetParameterAsText(4) ##Total building exposure
buildType = arcpy.GetParameterAsText(5) ##Building construction type - this influences loss ratios
outGDB = arcpy.GetParameterAsText(6)##Output Location (gdb)
buffName = arcpy.GetParameterAsText(7) ##Buffer name - this will persist for intersect name(s) too
BlastZones = arcpy.GetParameterAsText(8) ##Blast zones to use - DHS or Lloyds.
BlastSize = arcpy.GetParameterAsText(9) ##Blast size - drop down list
outCorName = arcpy.GetParameterAsText(10) ##output coordinate system - drop down list
lossTable = arcpy.GetParameterAsText(11) ##Output loss table (to GDB)
outShape = arcpy.GetParameterAsText(12) ##Output to shapefile - T or F
outCSV = arcpy.GetParameterAsText(13) ##Output loss table as csv - T or F

should be populating either a class, or maybe a dictionary. If you do it as a dictionary, then:

vars = {name: arcpy.GetParameterAsText(i)
 for i, name in enumerate((
 'Pfc', # Blast ID field - can be object ID or user defined/named field
 'idField', # ID of blast location to be used
 # ...
 ))
 }

Use OS methods

This seems harmless:

Pfc = str(desc.path) + "\\" + desc.name

though it could be cleaned up somewhat by the use of a raw f-string:

Pfc = rf'{desc.path}\{desc.name}'

But wait! Not all operating systems use a backslash as a path separator! So instead you should do

Pfc = os.path.join(desc.path, desc.name)

Just use literals

This:

TField1 = "DamageLevel"
TField2 = "ImpactZoneByFeet"
TField3 = "PropertyFootprint_Area"
TField4 = "PropertyImpacted_Area"
TField5 = "TotalExposure"
TField6 = "ImpactRatio"
TField7 = "ImpactExposure"
TField8 = "LossRatio"
TField9 = "Loss"
TablefList = [TField1,TField2,TField3,TField4,TField5,TField6,TField7,TField8,TField9]

should really just be

table_f_list = [
 'DamageLevel',
 'ImpactZoneByFeet',
 # ...
]

Parametrize hard-coded paths

If this never changes:

tableTemplate = r"O:\LossTable"

then it should be made capitalized (TABLE_TEMPLATE). However, it'd be better represented as an environmental variable or command-line argument.

Don't use comma notation for exceptions

This:

except Exception, e:

should be converted to the "new" (Python 2.6 and later, so... really not new at all) syntax:

except Exception as e:

Answer 3

Reading this code, two immediate things came to mind:

1. use multi-assignments. Instead of:

   for row in sCursor:
    bID = str(row[0])
    BR = int(row[1])
    intArea = row[2]
   

   try:

   for (bID, BR, intArea) in sCursor:
    bID, BR = str(bID), int(BR)
   

   and instead of:

   tmp = blastSize.rsplit(";")
   IZ1 = tmp[0]
   IZ2 = tmp[1]
   IZ3 = tmp[2]
   

   you can just do:

   IZ1, IZ2, IZ3 = blastSize.rsplit(";")
   

2. use dicts instead of chained if..then. Instead of:

   if BlastZones == "DHS":
   if BlastSize == "Backpack":
    blastSize = "26;80;125"
    LR1 = 1
    LR2 = 0.25
    LR3 = 0.1
   elif BlastSize == "truck1Ton":
    blastSize = "85;260;375"
    LR1 = 1
    LR2 = 0.25
    LR3 = 0.1
   elif BlastSize == "truck2Ton":
    blastSize = "100;335;500"
    LR1 = 1
    LR2 = 0.25
    LR3 = 0.1
   elif BlastSize == "truck5Ton":
    blastSize = "140;450;640"
    LR1 = 1
    LR2 = 0.25
    LR3 = 0.1
   elif BlastSize == "truck6Ton":
    blastSize = "175;590;980"
    LR1 = 1
    LR2 = 0.25
    LR3 = 0.1
   

   you can do:

   blastSizes = {"Backpack": "26;80;125",
    "truck1Ton": "85;260;375",
    "truck2Ton": "100;335;500",
    "truck5Ton": "140;450;640",
    "truck6Ton": "175;590;980"
    }
   if BlatZones == "DHS":
    LR1, LR2, LR3 = 1, 0.25, 1
    blastSize = blastSizes[BlastSize]
   

3. avoid empty statements. Instead of:

   if os.path.exists(outFolder):
    pass
   else:
    arcpy.CreateFolder_management(outGDB.rsplit("\\",1)[0],str(buffName+"_output"))
   

   use your boolean operators:

    if not os.path.exists(outFolder):
    arcpy.CreateFolder_management(outGDB.rsplit("\\",1)[0],str(buffName+"_output")) 
   

4. use introspection. Instead of:

   ##Set output coordinates
   if outCorName == "NAD83":
    outCor = NAD83
   elif outCorName == "WebMerc":
    outCor = WebMerc
   elif outCorName == "BNG":
    outCor = BNG
   

   you can do:

   outCor = locals()[outCorName]
   

• Reinderien: "use scatter/gather" That's not what you showed. Scatter/gather is the name of a vectored I/O strategy. What you're doing is tuple assignment.