Spatial analysis using ArcPy
Institute for Environmental and Spatial Analysis...University of North Georgia
Contents
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:
- The
Describe
function returns a Describe object with Layer properties. - The
da.Describe
function returns a Describe dictionary with Layer properties.
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 features
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
- the name of the urbanarea with the least hospitals found so far
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.3 Print selected fire stations for each high school
print(f'=== {school_name} ===')
with arcpy.da.SearchCursor(fire, ['LABEL']) as fire_cur:
for fire_row in fire_cur:
print(f' {fire_row[0]}')
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]}')