I've created a Python Script that determines the time to traverse a terrain based on Toblers Hiking Function. I initially created a sample set and clipped the DEM accordingly. The script ran the six school sites just over a minute. When I scaled up the process to include the full DEM, the same six school sites took just over an hour to process. I understand why this is happening as my Python script has to create a new cost raster for the entire study area for each school site.
What I looking for advice on is how I can I limit the processing extent for each school site based on the maximum walking time provided by the user. Based on Toblers Hiking Function the average walking distance over flat terrain is 5.037 km\hr. The idea behind Toblers Function is that as one traverses across a terrain the amount of distance one can cover decreases as the slope increases. Hence if the users maximum walking time is 60 min the maximum distance once can traverse if the terrain was flat would be 5.037 km. I'd like to set the extent to 5.037km using the site as the centroid of the processing extent. This would need to be updated for each site as to ensure the processing extent to be correctly set for each site.
I've attached Python Script below:
'''
Created on Jul 16, 2015
Calculate the time to traverse a terrain based on
Toblers Hiking Function based on slope
@author: PeterW
'''
# import system modules and site packages
import os
import arcpy
from arcpy.sa import *
import time
def hms_string(sec_elapsed):
h = int(sec_elapsed / (60 * 60))
m = int((sec_elapsed % (60 * 60)) / 60)
s = sec_elapsed % 60
return "{}h:{:>02}m:{:>05.2f}s".format(h, m, s)
# End hms_string
start_time = time.time()
# set environment settings
arcpy.env.overwriteOutput = True
# check out extensions
arcpy.CheckOutExtension("Spatial")
# set input and output workspace
inws = r"F:\Projects\2015\G111741\Model01\Rasters"
fgdb = r"F:\Projects\2015\G111741\Model01\Model01.gdb"
arcpy.env.workspace = inws
# calculate the slope in degrees
inras = "raw2"
out_measurement = "DEGREE"
# spatial reference for output results
coordsys = arcpy.Describe(inras).spatialReference
slope_deg = Slope(inras, out_measurement)
sites = os.path.join(fgdb,"Schools")
vertical_factor = arcpy.sa.VfTable(os.path.join(inws, "ToblerAway.txt")) #@UndefinedVariable
# create output feature class to store walk time polylines for each site
walk_poly = os.path.join(fgdb, "walk_poly")
# check if walk_poly exists else create it from scratch
if arcpy.Exists(walk_poly):
arcpy.Delete_management(walk_poly)
arcpy.CreateFeatureclass_management(fgdb, "walk_poly", "POLYLINE", "", "", "", coordsys, "", "", "", "")
arcpy.AddField_management(walk_poly, "Name", "TEXT", "", "", 25, "", "", "", "")
arcpy.AddField_management(walk_poly, "Contour", "SHORT", "", "", "", "", "", "", "")
else:
arcpy.CreateFeatureclass_management(fgdb, "walk_poly", "POLYLINE", "", "", "", coordsys, "", "", "", "")
arcpy.AddField_management(walk_poly, "Name", "TEXT", "", "", 25, "", "", "", "")
arcpy.AddField_management(walk_poly, "Contour", "SHORT", "", "", "", "", "", "", "")
# walking distance\time (1km = 12min; 2km = 24min; 2.5km = 30min; 5km = 60min)
walk_int = [12, 24, 30, 60]
with arcpy.da.SearchCursor(sites, ["Name"]) as cursor: #@UndefinedVariable
count = 0
for row in cursor:
count += 1
site_name = row[0].replace(" ","_")
sql_exp1 = "Name = '{}'".format(row[0])
site_lyr = site_name
arcpy.MakeFeatureLayer_management(sites, site_lyr, sql_exp1)
arcpy.env.extent = slope_deg
cost = PathDistance(site_lyr, slope_deg, "","","",inras, vertical_factor,"","")
print ("Processing {} Cost Raster").format(site_name)
costmin = Times(cost, 60)
print ("Processing {} Cost per Minute Raster").format(site_name)
walk_name = "walk_cont" + "_" + str(count)
walk_cont = os.path.join("in_memory", walk_name) # save to in_memory
print ("Processing & Sorting {} Walk Contours").format(site_name)
ContourList(costmin, walk_cont, walk_int)
walk_name2 = walk_name + "_" + "sorted"
walk_sorted = os.path.join("in_memory", walk_name2)
arcpy.Sort_management(walk_cont, walk_sorted, [["Contour", "ASCENDING"]],"")
arcpy.AddField_management(walk_sorted, "Name", "TEXT", "", "", 50, "", "", "", "")
sql_exp2 = "'{}'".format(row[0])
arcpy.CalculateField_management(walk_sorted, "Name", sql_exp2, "PYTHON_9.3")
with arcpy.da.SearchCursor(walk_sorted,["SHAPE@", "Name", "Contour"]) as scur: #@UndefinedVariable
with arcpy.da.InsertCursor(walk_poly, ["SHAPE@", "Name", "Contour"]) as icur: #@UndefinedVariable
for srow in scur:
icur.insertRow(srow)
print ("Completed Processing All Sites")
# check in extensions
arcpy.CheckInExtension("Spatial")
end_time = time.time()
print "It took {} to execute this".format(hms_string(end_time - start_time))
