Dr. Cho’s Website
Course Materials

Spatial analysis using ArcPy

Dr. Huidae Cho
Institute for Environmental and Spatial Analysis...University of North Georgia

1   Data set

In this lecture, we will perform various spatial analyses using ArcPy.

Add all the layers in the nc_spm_08_grass7_exercise data set.

2   How many hospitals are there in a selected urban area polygon?

2.1   Check if any features are selected

How do you know if there are selected features?

There are two ways:

Select one feature in urbanarea.

desc1 = arcpy.Describe('urbanarea')    # returns a Describe object
desc2 = arcpy.da.Describe('urbanarea') # returns a Describe dictionary

print(desc1.FIDSet)        # a string with the list of selected feature IDs separated by a space and semi-colon
print(desc2.get('FIDSet')) # a list of selected feature IDs

Select more features.

print(desc1.FIDSet)        # new set; dynamic
print(desc2.get('FIDSet')) # old set! snapshot; be careful!

print(arcpy.da.Describe('urbanarea').get('FIDSet')) # this is ok

fidset = [int(x) for x in desc1.FIDSet.split(';')] # create a list
print(fidset) # fidset is now similar to desc2.get('FIDSet'); both are static

They are different!

2.2   Allow only one feature selection and get that geometry

def get_selected_geometry(layer):
  desc = arcpy.Describe(layer)
  fidset = [int(x) for x in desc.FIDSet.split(';')] if desc.FIDSet else []
  if len(fidset) == 1:
    with arcpy.da.SearchCursor(layer, ['SHAPE@'], "{} = {}".format(desc.OIDFieldName, fidset[0])) as cur:
      geom = cur.next()[0]
  else:
    raise Exception('Select one feature!')
  return geom

2.3   Select hospitals within the selected urbanarea polygon

urban = 'urbanarea'
hosp = 'hospitals'

sel_urban_polygon = get_selected_geometry(urban)
arcpy.SelectLayerByLocation_management(hosp, 'WITHIN', sel_urban_polygon)

2.4   Count how many hospital features are selected

hosp_desc = arcpy.Describe(hosp)
print(len([int(x) for x in hosp_desc.FIDSet.split(';')] if hosp_desc.FIDSet else []))

2.5   Complete code

def get_selected_geometry(layer):
  desc = arcpy.Describe(layer)
  fidset = [int(x) for x in desc.FIDSet.split(';')] if desc.FIDSet else []
  if len(fidset) == 1:
    with arcpy.da.SearchCursor(layer, ['SHAPE@'], "{} = {}".format(desc.OIDFieldName, fidset[0])) as cur:
      geom = cur.next()[0]
  else:
    raise Exception('Select one feature!')
  return geom

urban = 'urbanarea'
hosp = 'hospitals'

# get the selected urban area polygon
sel_urban_polygon = get_selected_geometry(urban)

# select hospital features inside the polygon
arcpy.SelectLayerByLocation_management(hosp, 'WITHIN', sel_urban_polygon)

# see how many hospital features are selected
hosp_desc = arcpy.Describe(hosp)
print(len([int(x) for x in hosp_desc.FIDSet.split(';')] if hosp_desc.FIDSet else []))

3   Exercise: Find which urbanarea features have the least or no hospital features inside them

Hints

  • Loop through urbanarea features
  • Use two variables:
    • the name of the urbanarea with the least hospitals found so far
      • you may need a list because there can be multiple urbanarea features with the same least number of hospitals
    • the least number of hospitals

4   Find fire stations within a specified distance from all high schools

4.1   Search for all high schools

schools = 'schools_wake'
fire = 'firestations'
distance = 2000 # meters

with arcpy.da.SearchCursor(schools, ['SHAPE@', 'NAMELONG'], "GLEVEL = 'H'") as school_cur:

4.2   For each high school, create a buffer and select fire stations

  for school_row in school_cur:
    school_name = school_row [1]
    school_buf = school_row [0].buffer(distance)
    arcpy.SelectLayerByLocation_management(fire , 'WITHIN', school_buf)

4.4   Complete code

schools = 'schools_wake'
fire = 'firestations'
distance = 2000 # meters

with arcpy.da.SearchCursor(schools, ['SHAPE@', 'NAMELONG'], "GLEVEL = 'H'") as school_cur:
  for school_row in school_cur:
    school_name = school_row [1]
    school_buf = school_row [0].buffer(distance)

    arcpy.SelectLayerByLocation_management(fire , 'WITHIN', school_buf)

    print(f'=== {school_name} ===')

    with arcpy.da.SearchCursor(fire, ['LABEL']) as fire_cur:
      for fire_row in fire_cur:
        print(f'  {fire_row[0]}')