HiveBrain v1.2.0
Get Started
← Back to all entries
patternpythonMinor

ArcPy script to analyze land use in counties

Submitted by: @import:stackexchange-codereview··
0
Viewed 0 times
codereviewgeospatialstackoverflowfile-systemarcpypythonpython-2.x
scriptanalyzelandarcpycountiesuse

Problem

We have written a Python script using arcpy modules. It is was written by Python beginners and many parts of the code are written in 'unpythonic way'. The goal is to re-write or address 'unpythonic' code. The software used for this script are ArcMap 10.4 and written in Python 2.7. The complete script is available on GitHub.

The main directory is "C:/A__P6_GIS4/" and contains twelve sub-dirs the sub directories are actually names of counties, and each county sub-dir contains all relevant inputs inside the geodatabase or gdb file. First, the script loops and creates a .txt file with names of all the sub-dirs then uses this to loop through each sub-dir and execute set of functions for each county using set of if-else statements. The functions can be broadly broken by turf grass (TG) model, fractional (Frac) model, forest model and final model. The final model utilizes outputs from the aforementioned models and then creates the final model.

Only contains main() due to character limitation:

```
def main():

ALL_start_time = time.time()
#ALL_start_time = timeit.default_timer()

if arcpy.CheckExtension("Spatial") == "Available":
arcpy.AddMessage("Checking out Spatial")
arcpy.CheckOutExtension("Spatial")
else:
arcpy.AddError("Unable to get spatial analyst extension")
arcpy.AddMessage(arcpy.GetMessages(0))
sys.exit(0)

# Create List of Counties for Loop
arcpy.Delete_management("C:/GIS_temp/county_list.txt")
co_list = "C:/GIS_temp/county_list.txt"
county_list = open(co_list,"a") #Open the list
MainDIR = "C:/A__P6_GIS4/" #Directory where LAND USE geodatabases are located

for i in os.listdir(MainDIR):
i_base = os.path.basename(i)
county_list.write(i_base + "\n")
county_list.close()

# Start Looping Through County List
county_list = open(co_list,"r")
for i in county_list: #A text file is needed to ensure looping works!
CoName = i.strip("\n

Solution

Main

Use an helper function to separate unit of work. Namely, I would use a function performing the required operations on a county and call it from main.

def main(main_directory="C:/A__P6_GIS4/"):
    if arcpy.CheckExtension("Spatial") != "Available":
        arcpy.AddError("Unable to get spatial analyst extension")
        arcpy.AddMessage(arcpy.GetMessages(0))
        sys.exit(1)
    arcpy.AddMessage("Checking out Spatial")
    arcpy.CheckOutExtension("Spatial")

    for county in os.listdir(main_directory):
        manage_county(main_directory, county)


And that's all you need in your main. Separating this manage_county function into turf, forest, mixed and final subfunctions could be a good thing to do too.

A few things to note:

  • sys.exit(0) means there was no error, so better use an exit status of 1 to indicate an error;



  • os.listdir can be used directly to iterate over county directory, there in no need in using a file as buffer;



  • using a parameter with default value can help with reusability/maintenance as the function can easily be tested and such value is not burried within the code;



  • timing and debug printing can be delegated to helper functions/decorators, more on that later.



Manage_county

The main issue with the rest of the code, is the amount of redundant line of code one can read. Once again helper function can help reduce the amount of repetition. Loops are also a great way to perform the same operation on copious amount of filenames.

You should also take some time to remove useless variables, such as BEACH which is defined, printed, tested for existence, but nothing usefull is done with it.

You should also read PEP 8 and the official naming conventions to make the code look like Python code; and avoid needless abreviations in your variable names.

Strings management is also a mess: there is a lot of useless call to str as the variables it is applied to are already strings; os.path.join is mainly applied to a single string, thus it is just noise; and str.format should be prefered to string concatenation.

```
FEATURE_NAME_PATTERN = '{}/{}_{}'
FEATURE_1M_PATTERN = '{}/{}_{}_1m'

def create_directory(root, directory_name):
directory = os.path.join(root, directory_name)
if not arcpy.Exists(directory):
arcpy.CreateFolder_management(root, directory_name)
return directory

def create_geodatabase(root, filename):
file_name = os.path.join(root, filename)
created = False
if not arcpy.Exists(file_name):
arcpy.CreateFileGDB_management(root, filename)
created = True
return file_name

def manage_county(root, county_name):
county_directory = os.path.join(root, county_name)
output_directory = create_directory(county_directory, 'Outputs')
tiff_directory = create_directory(output_directory, county_name + '_FINAL')

imputs = os.path.join(county_directory, county_name + '_Inputs.gdb') # Former CoGDB
temp_1m, _ = create_geodatabase(output_directory, 'Temp_1m.gdb') # Former TempGDB
temp_10m, _ = create_geodatabase(output_directory, 'Temp_10m.gdb') # Former Temp10GDB
final_10m, _ = create_geodatabase(output_directory, 'Final_10m.gdb') # Former Final_10m
final_1m, created = create_geodatabase(output_directory, 'Final_1m.gdb') # Former LuGDB
if created:
for feature in ['IR', 'INR', 'TCoI', 'WAT']:
feature_in = FEATURE_1M_PATTERN.format(inputs, county_name, feature)
feature_out = FEATURE_1M_PATTERN.format(final_1m, county_name, feature)
arcpy.Copy_management(feature_in, feature_out)
arcpy.Copy_management(
FEATURE_NAME_PATTERN.format(inputs, county_name, 'LC'),
FEATURE_NAME_PATTERN.format(final_1m, county_name, 'LandCover'))

arcpy.env.overwriteOutput = True
coord_data = FEATURE_NAME_PATTERN.format(inputs, county_name, 'Snap')
ir_1m_path = FEATURE_1M_PATTERN.format(final_1m, county_name, 'IR')
arcpy.env.outputCoordinateSystem = arcpy.Describe(coord_data).spatialReference
arcpy.env.workspace = temp_1m
arcpy.env.scratchWorkspace = temp_1m
arcpy.env.extent = ir_1m_path
arcpy.env.parallelProcessingFactor = "100%"
arcpy.env.snapRaster = ir_1m_path #location of the default snap raster

#------------------------- TURF & FRACTIONAL MODELS -----------------------------
for parcel in ['IMP', 'IMP2']:
arcpy.Delete_management('{}/Parcel_{}'.format(temp_1m, parcel))
for feature in ['INRmask', 'RTmask', 'Parcels_TURFtemp', 'Parcels_TURF', 'TURF_parcels', 'Parcels_FTGtemp', 'Parcels_FTG', 'FTG_parcels', 'TGmask', 'FTGmask', 'TURFtemp', 'FTGtemp', 'FTGtemp2', 'FTGtemp3', 'FINRtemp']:
arcpy.Delete_management(FEATURE_NAME_PATTERN.format(temp_1m, county_name, feature))
for feature in ['TG', 'TCI', 'FTG1', 'FTG2', 'FTG3', 'FINR']:
arcpy.Delete_management(FEATURE_1M_PATTERN.format(final_1m, county_name, feature))

# Call each function, passing the nece

Code Snippets

def main(main_directory="C:/A__P6_GIS4/"):
    if arcpy.CheckExtension("Spatial") != "Available":
        arcpy.AddError("Unable to get spatial analyst extension")
        arcpy.AddMessage(arcpy.GetMessages(0))
        sys.exit(1)
    arcpy.AddMessage("Checking out Spatial")
    arcpy.CheckOutExtension("Spatial")

    for county in os.listdir(main_directory):
        manage_county(main_directory, county)
FEATURE_NAME_PATTERN = '{}/{}_{}'
FEATURE_1M_PATTERN = '{}/{}_{}_1m'


def create_directory(root, directory_name):
    directory = os.path.join(root, directory_name)
    if not arcpy.Exists(directory):
        arcpy.CreateFolder_management(root, directory_name)
    return directory


def create_geodatabase(root, filename):
    file_name = os.path.join(root, filename)
    created = False
    if not arcpy.Exists(file_name):
        arcpy.CreateFileGDB_management(root, filename)
        created = True
    return file_name


def manage_county(root, county_name):
    county_directory = os.path.join(root, county_name)
    output_directory = create_directory(county_directory, 'Outputs')
    tiff_directory = create_directory(output_directory, county_name + '_FINAL')

    imputs = os.path.join(county_directory, county_name + '_Inputs.gdb') # Former CoGDB
    temp_1m, _ = create_geodatabase(output_directory, 'Temp_1m.gdb') # Former TempGDB
    temp_10m, _ = create_geodatabase(output_directory, 'Temp_10m.gdb') # Former Temp10GDB
    final_10m, _ = create_geodatabase(output_directory, 'Final_10m.gdb') # Former Final_10m
    final_1m, created = create_geodatabase(output_directory, 'Final_1m.gdb') # Former LuGDB
    if created:
        for feature in ['IR', 'INR', 'TCoI', 'WAT']:
            feature_in = FEATURE_1M_PATTERN.format(inputs, county_name, feature)
            feature_out = FEATURE_1M_PATTERN.format(final_1m, county_name, feature)
            arcpy.Copy_management(feature_in, feature_out)
        arcpy.Copy_management(
            FEATURE_NAME_PATTERN.format(inputs, county_name, 'LC'),
            FEATURE_NAME_PATTERN.format(final_1m, county_name, 'LandCover'))

    arcpy.env.overwriteOutput = True
    coord_data = FEATURE_NAME_PATTERN.format(inputs, county_name, 'Snap')
    ir_1m_path = FEATURE_1M_PATTERN.format(final_1m, county_name, 'IR')
    arcpy.env.outputCoordinateSystem = arcpy.Describe(coord_data).spatialReference
    arcpy.env.workspace = temp_1m
    arcpy.env.scratchWorkspace = temp_1m
    arcpy.env.extent = ir_1m_path
    arcpy.env.parallelProcessingFactor = "100%"
    arcpy.env.snapRaster = ir_1m_path  #location of the default snap raster

    #------------------------- TURF & FRACTIONAL MODELS -----------------------------
    for parcel in ['IMP', 'IMP2']:
        arcpy.Delete_management('{}/Parcel_{}'.format(temp_1m, parcel))
    for feature in ['INRmask', 'RTmask', 'Parcels_TURFtemp', 'Parcels_TURF', 'TURF_parcels', 'Parcels_FTGtemp', 'Parcels_FTG', 'FTG_parcels', 'TGmask', 'FTGmask', 'TURFtemp', 'FTGtemp', 'FTGtemp2', 'FTGtemp3', 'FINRtemp']:
        arcpy.Delete_management(FEATURE_NAME_PATTERN.format(temp_1m, county_name, feature))
    for feature in ['TG', 'TCI', 'FTG1', 'FTG2', 'FTG3', 'FINR']:
        arcpy.Delete_management(FEATURE_1M_PATTERN.format(final_1m, county_name, feature))

    # Call each function, passing the necessary variables...
    turf_1(final_1m, county_name, temp_1m)
    turf_2(inputs, county_name, temp_1m)
 
def turf_1(final_1m, county_name, temp_1m):
    INRmask = FEATURE_NAME_PATTERN.format(temp_1m, county_name, 'INRmask')
    IR = FEATURE_1M_PATTERN.format(final_1m, county_name, 'IR')
    INR = FEATURE_1M_PATTERN.format(final_1m, county_name, 'INR')
    TCI = FEATURE_1M_PATTERN.format(final_1m, county_name, 'TCoI')
    ...
import time
from functools import wraps


def timer(func):
    @wraps(func)
    def wrapper(*args):
        start = time.time()  # or time.perf_counter() in Python 3
        print 'Starting', func.__name__, args
        func(*args)
        end = time.time()  # or time.perf_counter()
        print 'Computation time:', end - start
@timer
def manage_county(root, county_name):
    # rest of the code

Context

StackExchange Code Review Q#149591, answer score: 2

Revisions (0)

No revisions yet.