Performing a spatial join to obtain attribute field summaries of nested polygons that are located within polygon regions using ArcPy

I have created an ArcPy script tool that finds depressions from an input DEM, and then generates closed contours (both polyline and polygon outputs) as well as regions (raster and polygon outputs). I need my tool to be able to output a final layer that provides summary attributes of the closed contour polygons that fall within each region from the regions polygon feature class.

To do this, I know I need to perform a spatial join with the closed contour (cc) polygons and the regions polygons to obtain statistics from the attributes of the cc polygons within the regions; however, I want multiple statistics for each attribute field, such as ‘Max’,’Min’,and ‘Range’. I can successfully do this using the ‘Spatial Join’ within ArcMap 10.6’s user interface, but I want to do this using ArcPy – I think I need to use fieldmappings to do this and set up merge rules, but I’m not sure how to do this step. I have used fieldmappings to add both tables, and to remove fields that I don’t want to retain, but I am not sure where to proceed from there.

I have included a code snippet of my current fieldmappings before performing the spatial join below:

outFeatures = os.path.join(outLoc,outElim+ "_Final")
fieldmappings = arcpy.FieldMappings()

fieldmappings.addTable(outRegions)
fieldmappings.addTable(outElim)

keepers = ["area_sqft","perim_ft","Min_Z","Max_Z","Range_Z","circularityIndex","solidity","gridcode","length_ft","MBG_Width","MBG_Length","MBG_Orientation"]

#remove all output fields that aren't desired, then perform spatial joins
for field in fieldmappings.fields:
        if field.name not in keepers:
                fieldmappings.removeFieldMap(fieldmappings.findFieldMapIndex(field.name))

#set merge rules here to obtain min, max, and range for kept attribute fields


arcpy.SpatialJoin_analysis(outRegions,outElim,outFeatures,"JOIN_ONE_TO_ONE","KEEP_COMMON",fieldmappings,"CONTAINS","","")

#Set up fieldmappings for spatial join 
fieldNamesToMin = ["Min_Z","circularityIndex_min","solidity_min","MBG_Width_min","MBG_Length_min","MBG_Orientation_min"]
fieldNamesToMax = ["area_sqft","length_ft","Max_Z","circularityIndex_max","solidity_max","MBG_Width_max","MBG_Length_max","MBG_Orientation_max"]

fieldmappings = arcpy.FieldMappings()

fieldmappings.addTable(outRegions)
fieldmappings.addTable(outElimSummaryFields)

#remove undesired fields in final output
drops = ["perim_ft","gridcode"]

for field in fieldmappings.fields:
        if field.name in drops:
              
fieldmappings.removeFieldMap(fieldmappings.findFieldMapIndex(field.name))

#loop through field names for different merge rules (min and max values), and update field map with new merge rule, then perform spatial join
for fieldName in fieldNamesToMin:
        fieldIndexMin = fieldmappings.findFieldMapIndex(fieldName)
        fieldMap = fieldmappings.getFieldMap(fieldIndexMin)
        fieldMap.mergeRule = "MIN"

        fieldmappings.replaceFieldMap(fieldIndexMin,fieldMap)

for fieldName in fieldNamesToMax:
        fieldIndexMax = fieldmappings.findFieldMapIndex(fieldName)
        fieldMap2 = fieldmappings.getFieldMap(fieldIndexMax)
        fieldMap2.mergeRule = "MAX"

        fieldmappings.replaceFieldMap(fieldIndexMax,fieldMap2) 

outFeatures = os.path.join(outLoc,outElim+ "_Final")

arcpy.SpatialJoin_analysis(outRegions,outElimSummaryFields,outFeatures,"JOIN_ONE_TO_ONE","KEEP_COMMON",fieldmappings,"CONTAINS","","")