1

I am trying to spatially join characteristics from a larger unit (here: county) to a smaller unit (say, local district), such that for every local district I will also know in which county it is in.

The problem arises when a district lies in two (or more) different counties. In this case, I would like the local district to take the attributes from the county that covers most of the district's area. In the screenshot below,

enter image description here

the two local districts marked yellow were assigned attributes from the Greene County despite the fact that most of their area is within the Sumter County. The reason is that I used standard code

target_features = "district.shp"
join_features = "county.shp"
arcpy.SpatialJoin_analysis(target_features, join_features, out_feature_class)

and without any further options specified ArcPy (I guess) simply used the first county information it got hold off.

Is there an option that allows me to account for the area that a district shares with a county?

I have looked at the match_option CLOSEST but I think it does not do the job, because - depending on Green County's total size - its center might very well be closer to a borderline district even if that should be part of Sumter County (like the two yellow districts above).

I am new to ArcPy/GIS.

3
  • 1
    A spatial join, by design, does not fracture features on boundaries. You need to look in the Overlay toolset fof Union or Intersect
    – Vince
    Commented Nov 17, 2016 at 11:40
  • Thank you, Vince, good suggestion. I will look into union and intersect and see if they can help. Commented Nov 17, 2016 at 14:53
  • gis.stackexchange.com/questions/217729/…
    – FelixIP
    Commented Nov 17, 2016 at 18:41

2 Answers 2

1

You'll want to use intersect, and then find the output feature with the largest area per district.

#distrcit feature class
distFc = r"Path\to\Districts"
#district name field
distFld = "DISTRICT"

#county feature class
countyFc = r"Path\to\Counties"
#county name field
countyFld = "COUNTY"

#-----------

print "import"
import arcpy

print "intersecting districts and counties"
intersectFc = "C\Path\to\intersect\output"
arcpy.Intersect_analysis ([distFc, countyFc], intersectFc)

#size dictionary
sizeDi = {}

#county dictionary
countyDi = {}

print "iterating intersect feature class"
with arcpy.da.SearchCursor (intersectFc, ["SHAPE@AREA", distFld, countyFld]) as curs:
    for area, dist, county in curs:
        #try to get the area of the last intersect feature
        try:
            lastArea = sizeDi [dist]
        #if none, set last area to zero
        except KeyError:
            lastArea = 0
        #check if current area is greater than last area
        if area > lastArea:
            #current area is greater than last area
            #set district to county
            countyDi [dist] = county
            #set new size in size di
            sizeDi [dist] = area

"""
'countyDi' now has county the district falls into the most
key: district name
value: highest spatial overlap county name
"""
0

As commented by @Vince:

A spatial join, by design, does not fracture features on boundaries. You need to look in the Overlay toolset for Union or Intersect

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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