1

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?

enter image description here

3
  • 1
    Lots of possible issues here, including geometry errors. You should always reflexively Check Geometry whenever you get unexpected overlay results.
    – Vince
    Jun 10 at 12:09
  • Thanks for your answer! But what do you mean with checking the geometry? I think the geometry is pretty clear here. It's a single part feature and the boundaies are selected. The Cordinate Systems are also consistent. I used the method "centroid within" instead of "largest overlap" and this works pretty good.
    – Radde1683
    Jun 10 at 12:33
  • Check Geometry is a tool that looks for topology errors in geodata (Repair Geometry fixes them). Those errors cause tools that assume correct topology to fail.
    – Vince
    Jun 10 at 13:18

1 Answer 1

2

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
1
  • Thank you really much! This worked way better than the original tool.
    – Radde1683
    Jun 16 at 9:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.