I have a Polygon Shapefile that represents US counties and a Point Shapefile that represents hospitals. The hospital attributes table includes fields like "Number of Beds", "Number of Patients" and "Number of Doctors". I want to find the data for the hospital which has the maximum number of doctors in each county.

Lets say these Points fall inside county xyz.


1 ................ 20 ...... 200 .......... 10

2 ................ 12 ...... 300 .......... 5

3 ................ 15 ...... 100 .......... 13

My desired result should look like:


xyz ............3 ................ 15 ...... 100 .......... 13

Getting just maximum of Doctors is not what I need (or maximum of all fields).


2 Answers 2


Note: Edited answer based on additional info provided by OP

Try the following:

1) Combine the county and hospitals datasets, either via a spatial join or with identify (I personally prefer identify so I have a new layer to work with, but either should work). This should give you a single point dataset with all your hospital info, but with the addition of the county name that it lies in.

2) Run the SORT tool to create a sorted dataset based on the county name field and the Doctor field, each sorted in decending order. If you are likely to have multiple hospitals in a county with the same number of doctors, and you have another field which would be the deciding factor in determining which hospital wins in event of that tie (ex: if 2 hospitals have the same max number of doctors, chose the one of those with the greater number of beds), then add that field as well.

3) Run the Summary Statistics tool. For the STATISTICS FIELD(S) add each field from the hospital data you want data from (ID, doctors, patients, beds), with Statistics Type = First. For CASE FIELD, select the county name (or ID). This will esentially give you the first feature for each county, and since you already sorted it to put the hospital with the max doctors at the top for each county, this should get you what you need.

  • Intersect would also work - Identity requires an Advanced License. The SORT option is a creative approach, particularly to deal with ties, but Summary Statistics offers a MAX statistic. You could leave out step 2 and just run SS on the table with county ID as case field and doctors as statistic with MAX as the type. In theory you'll know if you have a tie anywhere because the COUNT field will be greater than one.
    – Chris W
    Oct 24, 2014 at 20:28
  • That was basically my original response, but onfortunately that doesn't retain the other informaition for the hospital that has the max doctors. And thanks for the note about Identidy requiring an Advanced Licence - sometimes I forget that some of the tools I can use aren't available to everyone!
    – Beck
    Oct 24, 2014 at 20:38
  • True, but I suppose it depends on the desired output. SS generates a table, not geometry, so if you want to identify those hospitals and then select them or export them to a new dataset or anything, you'll have to join the result back using hospital id or something anyway. If the goal is just the tabular data, you're right, it has to be worked into the SS tool in some way.
    – Chris W
    Oct 24, 2014 at 21:45

Here's an arcpy solution. Run it within the ArcMap Python window. You need to change the three variables at the top (counties, hospitals, and field) to match your data.

The output is an in-memory copy of the counties layer called Counties_output with all the fields from the hospital layer joined to it for the maximum value. If there are multiple occurrences of the maximum value, there will be duplicate features for each one.

import arcpy

#Set these as the name of the layer in ArcMap's TOC
counties = "County_layer"
hospitals = "Hospital_layer"
#Set this as the name of the field to find max value for
field = "Doctors"

object_id = []
hospital_fl = arcpy.MakeFeatureLayer_management(hospitals, "hospital_fl")

with arcpy.da.SearchCursor(counties, ("SHAPE@")) as rows:
    for row in rows:
        polygon_geom = row[0]
        arcpy.SelectLayerByLocation_management(hospital_fl, "COMPLETELY_WITHIN", polygon_geom)
        doctors = {row[0]: row[1] for row in arcpy.da.SearchCursor(hospital_fl, ("OID@", field))}
        max_doctors = max(doctors.values())
        for k, v in doctors.items():
            if v == max_doctors:

arcpy.SelectLayerByAttribute_management(hospital_fl, "NEW_SELECTION", "OBJECTID IN {0}".format(tuple(object_id)))
arcpy.SpatialJoin_analysis(counties, hospital_fl, "in_memory\\Counties_output", "JOIN_ONE_TO_MANY", "KEEP_ALL", "#", "CONTAINS")

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.