Creating a Multi-Step Spatial Join Tool Script (One point file to many polygons) with one point output

Question

I am wanting to create a ArcPy Tool Script that will join a point file to multiple polygons. Ideally, the point file (input_fc) would be a parameter so it can be interchangeable. Then the polygons would be in one geodatabase (CensusGDB)(the files will be Census Blocks, Block Groups, Tracts and Counties). The output geodatabase would also be a parameter (out_gdb).

The final output would be one point file with just the key field for each of these files that were joined. This is what I have so far, doing the first spatial join, but will need the output of this first spatial join output to be the input for the next spatial joins.

##Spatial Join Census geographies
import arcpy
import os
CensusGDB = r"C:\Census_All.gdb"
out_gdb = arcpy.GetParameterAsText(1)
#Input the point feature class
input_fc = arcpy.GetParameterAsText(0)
out_fc = os.path.join(out_gdb, 'name')
#List the featureclasses within the Census gdb
arcpy.env.workspace = CensusGDB
featureclasses = arcpy.ListFeatureClasses()
# for fc in featureclasses:
for fc in featureclasses:
desc = arcpy.Describe(fc)
fc_path =os.path.join(CensusGDB,desc.name)
out_fc = os.path.join(out_gdb,desc.name)
    ##Start Spatial Joins
arcpy.analysis.SpatialJoin(target_features =input_fc,join_features=fc_path,output_feature_class = out_fc)

Any ideas on if this is the right way to go or if I need to try something else?

Answer 1

Assuming I understand your question, this worked for me with some mock data. Basically, I merged all the input data together and did a single Spatial Join, retaining only the specified fields.

I stored the names of my polygon feature classes and the associated "key" field that should be joined to the points as a Python dictionary. Not sure what the best way to take that input would be for an ArcGIS toolbox (perhaps a ValueTable object), as I typically work from an IDE without constructing the full toolboxes.

#IMPORTS----------------------
import arcpy as ap

import sys
import traceback

#SETTINGS----------------------
ap.env.workspace = r'C:\Users\vce50\Documents\ArcGIS\Projects\SCRATCH_Point_Join\SCRATCH_Point_Join.gdb'
ap.env.overwriteOutput = True

#INPUTS----------------------
in_pnt_fc = 'Point'

# Dictionary of feature classes to join, and their "key" field.
joining_fcs = {'County_Boundaries': 'COUNTY_FIP',
               'Municipal_Boundaries': 'MUNICIPALI',
               'Zip_Codes': 'ZIPCODE1'}

#MAIN----------------------
try:
    # Create Field Mappings of only the required fields.
    fieldMappings = ap.FieldMappings()
    for fc, fld in joining_fcs.items():
        fieldMap = ap.FieldMap()
        fieldMap.addInputField(fc, fld)
        fieldMappings.addFieldMap(fieldMap)
    print('Field Mappings built.')

    # Merge all the input datasets. Store output in memory.
    merged_fcs = ap.management.Merge(inputs=list(joining_fcs.keys()),
                        output=r'memory\Merged_FCs',
                        field_mappings=fieldMappings)
    print('Inputs merged.')

    # Spatial Join attribute values from merged FC to new Point FC.
    pnt_join = ap.analysis.SpatialJoin(target_features=in_pnt_fc,
            join_features=merged_fcs,
            out_feature_class=f'{in_pnt_fc}_Joined',
            field_mapping=fieldMappings)
    print('Polygon attributes spatially joined to point FC.')
    
    # Delete some fields added during Spatial Join.
    ap.management.DeleteField(pnt_join, ['Join_Count', 'TARGET_FID'])
    print('Fields dropped.')

except:
    etype, evalue, tback = sys.exc_info()
    tback_info = traceback.format_tb(tback)[0]
    err_msg = (f"Traceback Info:\n{tback_info}\n{etype.__name__}: {evalue}")
    print(err_msg)

My output points now contain the county ID, municipality, and zip.