Using Arcpy with multiprocessing – Part 3

The following post builds upon the script and methods developed in Part 1 and Part 2, so read them first!

The purpose of this series is not to give you a one line example that you can copy and paste to your code, but step through the process and make the underlying principles clear. Please note that Esri also has a blog post describing use of Python’s Multiprocessing library with Arcpy (in particular Spatial Analyst) which did not work for me in a complex situation involving Network Analyst, but is worth checking out if you are doing intensive tasks with Spatial Analysis.

The following Arcpy scripts refer to two GDBs which I have zipped into an archive; click here to download the archive (opens in new tab). If using these scripts from the command line or a Python IDE remember to change the paths to the Feature Classes!

Now that the most calculation intensive part of the process has been made into a parallelisable function we can actually implement the multiprocessing bit. Here I will describe basic use cases for two multiprocessing libraries; the inbuilt Python Multiprocessing, and Parallel Python. I include both as neither seems to work for all situations; Parallel Python has caused problems when ArcGIS extensions are checked out and I cannot get Multiprocessing to work when geocoding addresses. If you have odd problems with one library I would suggest trying the other.

Both libraries can use a Python list as the data structure to which jobs are appended, and from which the result is extracted. Each job is processed in turn by the first available worker (CPU/core). So, if there are six jobs in the list on a system with two CPUs, the first two jobs will start near simultaneously, then as each CPU finishes its task and becomes free it will start working on the next job in the list. This is known as asynchronous processing, as the order and timing of completion are unimportant. Both methods can use all available CPUs on the system, or they can be limited to a predefined number (see below), or they may be limited by the method itself – in this example, as there are only 6 Polygon types, a maximum of only 6 CPUs would be used.

'''
Step 4 of 4: Example of using Multiprocessing/Parallel Python with Arcpy...

Can be run either:
 1. from the command line/a Python IDE (adjust paths to feature classes, as necessary)
 2. as a Script tool within ArcGIS (ensure 'Run Ptyhon script in Process' is NOT checked when importing)

 The Parallel Python library must be installed before it can be used.
'''

import arcpy
import multiprocessing
import time
try:
	import pp
	forceMP = False
except ImportError:
	forceMP = True

def performCalculation(points_fC, polygons_fC, searchDist, typeList, calcPlatform_input=None):
	'''performCalculation(pointsFeatureClass, polygonsFeatureClass, searchDistance, typeList, calcPlatform)

	All inputs are specific.

	calcPlatform is either 'mp', to use the inbuilt Multiprocessing library, or 'pp', to use Parallel Python (if it exists, otherwise defaults to mp)

	'''

	# ---------------------------------------------------------------------------
	## Pre-calculation checks
	# ---------------------------------------------------------------------------

	# Check calculation platform is valid, or resort to default...
	defaultCalcTpye = 'mp'
	calcTypeExplainDict = {'pp':'Parallel Python', 'mp':'Multiprocessing'}
	if forceMP: # unable to import Parallel Python library (it has to be installed seperately)
		arcpy.AddMessage("   WARNING: Cannot load Parallel Python library, forcing Multiprocessing library...")
		calcPlatform = 'mp'
	elif (calcPlatform_input not in ['mp', 'pp']) or (calcPlatform_input == None): # Invalid/no input, resort to default
		arcpy.AddMessage("   WARNING: Input calculation platform '%s' invalid; should be either 'mp' or 'pp'. Defaulting to %s..." % (calcPlatform_input, calcTypeExplainDict[defaultCalcTpye]))
		calcPlatform = defaultCalcTpye
	else:
		calcPlatform = calcPlatform_input

	# ---------------------------------------------------------------------------
	## Data extraction (parallel execution)
	# ---------------------------------------------------------------------------

	searchDist = int(searchDist) # convert search distance to integer...

	# check all datasets are OK; if not, provide some useful output and terminate
	valid_points = arcpy.Exists(points_fC)
	arcpy.AddMessage("Points Feature Class: "+points_fC)
	valid_polygons = arcpy.Exists(polygons_fC)
	arcpy.AddMessage("Polygons Feature Class: "+polygons_fC)
	dataCheck  = valid_points & valid_polygons

	# Check if Inputs are valid and print message if bad, else start calculation...
	if not dataCheck:
		if not valid_points:
			arcpy.AddError("Points database or feature class, %s,  is invalid..." % (points_fC))
		if not valid_polygons:
			arcpy.AddError("Polygons database or feature class, %s,  is invalid..." % (polygons_fC))

	else:
		for type in typeList: # add fields to Points file
			arcpy.AddField_management(points_fC, "polygon_type%s_Sum" % type, "DOUBLE")
			arcpy.CalculateField_management(points_fC, "polygon_type%s_Sum" % type, 0, "PYTHON")
			arcpy.AddField_management(points_fC, "polygon_type%s_Count" % type, "DOUBLE")
			arcpy.CalculateField_management(points_fC, "polygon_type%s_Count" % type, 0, "PYTHON")
			arcpy.AddMessage("    Added polygon_type%s_Sum and polygon_type%s_Count fields to Points." % (type,type))

		# dictionary: pointsDataDict[pointID][Type]=[sum of Polygon type weightings, count of Ploygons of type]
		pointsDataDict = {}
		# jobs are added to this list, then processed by the available workers (CPUs)
		jobs = [] 

		if calcPlatform == 'mp':
			arcpy.AddMessage("    Utilising Python Multiprocessing library:")
			# initate the processing pool - use all available resources
			pool = multiprocessing.Pool()
			## pool = multiprocessing.Pool(2) # Example: limit the processing pool to 2 CPUs...
			for type in typeList:
				arcpy.AddMessage("      Passing %s to processing pool..." % type)
				# send jobs to be processed
				jobs.append(pool.apply_async(findPolygons, (points_fC, polygons_fC, type, searchDist))) 

			# collect results from the job server (waits until all the processing is complete)
			for job in jobs:
				# first output becomes the new dictionary
				if len(pointsDataDict.keys()) == 0:
					pointsDataDict.update(job.get())
				else:
					# later outputs should be added per key...
					for key in job.get():
						pointsDataDict[key].update(job.get()[key])
			del jobs

		elif calcPlatform == 'pp':
			ppservers=()
			# initate the job server - use all available resources
			job_server = pp.Server(ppservers=ppservers)
			## job_server = pp.Server(2, ppservers=ppservers) # Example: limit the processing pool to 2 CPUs...
			arcpy.AddMessage("    Utilising Parallel Python library:")
			for type in typeList: # For this specific example, it would only initate three processes anyway...
				arcpy.AddMessage("      Passing %s to processing pool..." % type)
				# send jobs to be processed
				jobs.append(job_server.submit(findPolygons, (points_fC, polygons_fC, type, searchDist), (), ("arcpy",))) 

			# collect results from the job server (waits until all the processing is complete)
			for job in jobs:
				if len(pointsDataDict.keys()) == 0: # first output becomes the new dictionary
					pointsDataDict.update(job())
				else:
					for key in job():
						pointsDataDict[key].update(job()[key]) # later outputs should be added per key...
			del jobs

		# ---------------------------------------------------------------------------
		## Writing data back to points file
		# ---------------------------------------------------------------------------

		arcpy.AddMessage("    Values extracted; writing results to Points...")
		pointsRowList = arcpy.UpdateCursor(points_fC)
		for pointsRow in pointsRowList:
			pointID = pointsRow.getValue("PointID")
			for type in typeList:
				pointsRow.setValue("polygon_type%s_Sum" % type, pointsRow.getValue("polygon_type%s_Sum" % type) + pointsDataDict[pointID][type][0])
				pointsRow.setValue("polygon_type%s_Count" % type, pointsRow.getValue("polygon_type%s_Count" % type) + pointsDataDict[pointID][type][1])

			pointsRowList.updateRow(pointsRow)
			del pointsRow

		del pointsRowList
		# just make sure any locks are cleared...
		del calcPlatform, calcPlatform_input, calcTypeExplainDict, dataCheck, defaultCalcTpye, job, key, pointID, pointsDataDict, points_fC, polygons_fC, type, valid_points, valid_polygons, searchDist, typeList

def findPolygons(points_FC, polygons_FC, Type, search_dist):
	funcTempDict = {}
	arcpy.MakeFeatureLayer_management(polygons_FC, '%s_%s' % (polygons_FC, Type))
	arcpy.MakeFeatureLayer_management(points_FC, '%s_%s' % (points_FC, Type))
	pointsRowList = arcpy.SearchCursor('%s_%s' % (points_FC, Type))
	for pointRow in pointsRowList: # for every origin
		pointID = pointRow.getValue("PointID")

		try:
			funcTempDict[pointID][Type] = [0,0]
		except KeyError: # first time within row
			funcTempDict[pointID] = {}
			funcTempDict[pointID][Type] = [0,0]

		arcpy.SelectLayerByAttribute_management('%s_%s' % (points_FC, Type), 'NEW_SELECTION', '"PointID" = \'%s\'' % pointID)
		arcpy.SelectLayerByLocation_management('%s_%s' % (polygons_FC, Type), 'INTERSECT', '%s_%s' % (points_FC, Type), search_dist, 'NEW_SELECTION')
		arcpy.SelectLayerByAttribute_management('%s_%s' % (polygons_FC, Type), 'SUBSET_SELECTION', '"polyType" = %s' % Type)
		polygonsRowList = arcpy.SearchCursor('%s_%s' % (polygons_FC, Type))
		for polygonsRow in polygonsRowList:
			funcTempDict[pointID][Type][0] += polygonsRow.getValue("polyWeighting")
			funcTempDict[pointID][Type][1] += 1

	return funcTempDict

if __name__ == '__main__':
	# Read the parameter values:
	#  1: points feature class (string to location, e.g. c:\GIS\Data\points.gdb\points01)
	#  2: polygons feature class (string to location)
	#  3: search distance for polygons, integer, e.g 500
	#  4: calculation type ('mp' to use Multiprocessing library, 'pp' to use Parallel Python library (if available, otherwise defaults to mp))
	pointsFC = arcpy.GetParameterAsText(0) # required
	polygonsFC = arcpy.GetParameterAsText(1) # required
	search_Distance = arcpy.GetParameterAsText(2) # required
	calcType = arcpy.GetParameterAsText(3) # optional (will default to MP)

	t_start = time.clock()

	type_list = [1,2,3,4,5,6] # polygon types to search for...

	# run calculation
	if '' in [pointsFC, polygonsFC, search_Distance]:# if not running from Arc, the input parameters all come out as ''
		pointsFC = "c:\\Work\\GIS\\Data\\Points.gdb\\points01"
		polygonsFC = "c:\\Work\\GIS\\Data\\Polygons.gdb\\polygons01"
		search_Distance = 1000
		performCalculation(pointsFC, polygonsFC, search_Distance, type_list)
	else:
		performCalculation(pointsFC, polygonsFC, search_Distance, type_list, calcType)

	arcpy.AddMessage("      ... complete in %s seconds." % (time.clock() - t_start))

When importing to use it as an Arc script tool, ‘Run Python script in Process’ must be unchecked. If you don’t do this, the script will try to open new instances of ArcMap/ArcCatalog and fail.

Execution times:
On my computer, the scripts posted in Part 2 all take about 300 seconds to execute if run from the command line, and about 400 seconds if run as a tool within ArcMap or ArcCatalog. The script above takes about 170 seconds to execute if run from the command line, and 340 seconds if run as a tool within ArcMap or ArcCatalog. The first set of processes through the queue (i.e. two if two CPUs are being used, four if four CPUs are being used) are slow, but later iterations run a lot faster. This is called overhead, and can often mean that parallelised operations will not run as fast as expected (i.e. going from one process to four processes does not speed up the total execution time by a factor of four). The overall effect of overhead obviously reduces as the number of iterations increases.

Feel free to ask questions or point out improvements with regard to any of these posts. The following is a list of tips/issues I have noticed with regards to Arcpy scripting, and multiprocessed applications in particular (note: multiporcessing refers to the general concept, Multiprocessing refers to the inbuilt Python library).

Notes on using multiprocessing:

1. Creating and debugging complex Arcpy scripts is a real pain in the ass, and it is even worse when you are using multiprocesssing. I wouldn’t recommend doing it unless it is really necessary and you will be able to use the code extensively in the future. Errors are a little harder to debug, as error message output isn’t produced until the job result is being accessed.
2. Python Multiprocessing vs. Parallel Python:There doesn’t seem to be much speed difference between the two, however:
   1. Parallel Python requires an extra package to be installed (Multiprocessing is included with Python)
   2. Parallel Python will fail if you use arcpy.addMessage() or arcpy.addError() within the parallelised function (but you can just pass back a string and print it from the result)
   3. Parallel Python has issues if you are using tools that must be ‘Checked Out’ before execution (such as Network Analyst)
   4. I have not been able to get the Multiprocessing library to work with Geocoding
3. arcpy.Exists() is practically essential – it should be used to check every input dataset. Not just because it will catch input errors, but without it, Arc will regularly fail to do basic operations for no apparent reason (I assume this is related to data locks). It may be that checking for existence forces a refresh that clears locks, but it is certainly not foolproof either.
4. Try…Except– I have found that using these statements, as ESRI seem to like doing in the Arcpy help/tutorials where they wrap the entire script in a block, will often cause the script to hang indefinitely when running from Arc. The best thing to do is test rigorously in Python (passing the paths as strings and printing them), where you get useful outputs, then try the script out in Arc. If it ran fine directly, but has problems at this stage, make sure to check that the paths coming out of Arc are what you were actually expecting by using print.
5. Cancelling a multiprocessing script running in Arc does not kill all the Python processes – unlike hitting Ctrl-C when running from the command line. Use task manager to find them all and end them.

If you are doing Python development, you may be interested in my Windows Dev Stack, which describes my development environment from high level technologies down to specific apps, and how they all work together.