5
\$\begingroup\$

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
\$\endgroup\$
1
  • \$\begingroup\$ 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 \$\endgroup\$ Commented Jan 7, 2019 at 12:38

3 Answers 3

Reset to default
1
\$\begingroup\$

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)
\$\endgroup\$
2
\$\begingroup\$

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:
\$\endgroup\$
1
\$\begingroup\$

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]
\$\endgroup\$
1
  • \$\begingroup\$ "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. \$\endgroup\$ Commented Jan 6, 2019 at 0:32

You must log in to answer this question.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.