Spatial Join with largest overlap does not work in ArcGIS Pro

Question

I'm trying to do a spatial join with a parcel shapefile and a soil-water shapefile. I tried to do a one-to-one spatial join with the match option "largest overlap", so the parcel shapefile gets the soil-water values of the areas it has the largest overlap with. But it only really works in areas where it has a 100% overlap.

I hope you can see it on the picture. The object in the red circle only has a 10% overlap with the value "mittel" from the water shapefile and a 90% overlap with the value "hoch". Still the shapefile that is created gets the value "mittel".

Did anyone else ever had the same problem and know a solution for that?

Answer 1

This toolbox worked well for us:

https://www.arcgis.com/home/item.html?id=e9cccd343bf84916bda1910c31e5eab2

In case the link breaks in the future, here's a .pyt with the same content:

import arcpy
import os
from arcpy import env
from arcpy.sa import *

class Toolbox(object):
    def __init__(self):
        self.label = "Largest Overlap Toolbox"
        self.alias = ""
        self.tools = [SimpleSpatialJoinLargestOverlap]

class SimpleSpatialJoinLargestOverlap(object):
    def __init__(self):
        self.label = "Simple Spatial Join (Largest Overlap)"
        self.description = "Creates a spatial join with the features that have the largest overlap"
        self.canRunInBackground = False

    def getParameterInfo(self):
        #Define parameter definitions
           
        out_work = arcpy.Parameter(
            displayName="Output Workspace",
            name="out_work",
            datatype="DEWorkspace",
            parameterType="Required",
            direction="Input")
        
        target_features = arcpy.Parameter(  #0
            displayName="Target Features",
            name="target_features",
            datatype="GPFeatureLayer",
            parameterType="Required",
            direction="Input")
        

        join_features = arcpy.Parameter(  #1
            displayName="Join Features",
            name="join_features",
            datatype="GPFeatureLayer",
            parameterType="Required",
            direction="Input")
        
        
        out_fc = arcpy.Parameter(  #2
            displayName="Output Feature Class",
            name="out_fc",
            datatype="GPFeatureLayer",
            parameterType="Required",
            direction="Output")
            
            
        keep_all = arcpy.Parameter( # 3
            displayName="Keep All",
            name="keep_all",
            datatype="GPBoolean",
            parameterType="Optional",
            direction="Input")
        keep_all.value = True   
        
       
        spatial_rel= arcpy.Parameter( # 4
            displayName = "largest_overlap",
            name = "spatial_rel",
            datatype = "GPString",
            parameterType = "Optional",
            direction = "Input")
        spatial_rel.filter.type = "ValueList"
        spatial_rel.filter.list = ['largest_overlap'] 
    
        params =   [target_features, #0
                    join_features,   #1
                    out_fc,          #2
                    keep_all,        #3
                    spatial_rel]     #4
                    
        return params            

    def isLicensed(self):
        return True

    def execute(self, parameters, messages):
        arcpy.env.overwriteOutput = True

        #Define parameters    
        target_features = parameters[0].valueAsText
        join_features = parameters[1].valueAsText
        out_fc = parameters[2].valueAsText
        keep_all = parameters[3].valueAsText
        spatial_rel = parameters[4].valueAsText
    
    
        if spatial_rel == "largest_overlap":
            # Calculate intersection between Target Feature and Join Features
            intersect = arcpy.analysis.Intersect([target_features, join_features], "in_memory/intersect", "ONLY_FID")
            # Find which Join Feature has the largest overlap with each Target Feature
            # Need to know the Target Features shape type, to know to read the SHAPE_AREA oR SHAPE_LENGTH property
            geom = "AREA" if arcpy.Describe(target_features).shapeType.lower() == "polygon" and arcpy.Describe(join_features).shapeType.lower() == "polygon" else "LENGTH"
            fields = ["FID_{0}".format(os.path.splitext(os.path.basename(target_features))[0]),
                    "FID_{0}".format(os.path.splitext(os.path.basename(join_features))[0]),
                    "SHAPE@{0}".format(geom)]
            overlap_dict = {}
            with arcpy.da.SearchCursor(intersect, fields) as scur:
                for row in scur:
                    try:
                        if row[2] > overlap_dict[row[0]][1]:
                            overlap_dict[row[0]] = [row[1], row[2]]
                    except:
                        overlap_dict[row[0]] = [row[1], row[2]]
    
            # Copy the target features and write the largest overlap join feature ID to each record
            # Set up all fields from the target features + ORIG_FID
            fieldmappings = arcpy.FieldMappings()
            fieldmappings.addTable(target_features)
            fieldmap = arcpy.FieldMap()
            fieldmap.addInputField(target_features, arcpy.Describe(target_features).OIDFieldName)
            fld = fieldmap.outputField
            fld.type, fld.name, fld.aliasName = "LONG", "ORIG_FID", "ORIG_FID"
            fieldmap.outputField = fld
            fieldmappings.addFieldMap(fieldmap)
            # Perform the copy
            arcpy.conversion.FeatureClassToFeatureClass(target_features, os.path.dirname(out_fc), os.path.basename(out_fc), "", fieldmappings)
            # Add a new field JOIN_FID to contain the fid of the join feature with the largest overlap
            arcpy.management.AddField(out_fc, "JOIN_FID", "LONG")
            # Calculate the JOIN_FID field
            with arcpy.da.UpdateCursor(out_fc, ["ORIG_FID", "JOIN_FID"]) as ucur:
                for row in ucur:
                    try:
                        row[1] = overlap_dict[row[0]][0]
                        ucur.updateRow(row)
                    except:
                        if not keep_all:
                            ucur.deleteRow()
            # Join all attributes from the join features to the output
            joinfields = [x.name for x in arcpy.ListFields(join_features) if not x.required]
            arcpy.management.JoinField(out_fc, "JOIN_FID", join_features, arcpy.Describe(join_features).OIDFieldName, joinfields)        
        
        
        return