I have a points shapefile and a polygon shapefile, each with millions of features. All points lie spatially within polygons so that there is one or more points within each polygon.

The points have an attribute called DOY. For each polygon, I want to find and extract the point with the minimum DOY to a new layer. So, the resulting points layer would have one point per polygon.

So far I have used Spatial Join to join the two shapefile attributes, but I am not sure how to extract the points I want. I was able to create a layer with a new field that tells me the minimum DOY value for each polygon, but I want to actually extract these points as well so I can identify them. What approach or geoprocessing tools do you recommend for this in ArcGIS?


1 Answer 1


Intersect your points and polygons. Then:

Option 1: If you have Advanced license one way is to Sort by point value then Delete Identical by polygon id field.

Option 2: Use arcpy:

import arcpy

#Adjust these three lines to match you input layer and fields
fc = r'points_inter_grid' #The intersection of your points and polygons output layer
polygon_id_field = 'ruta' #Polygon identifier field
point_value_field = 'DOY' #Point value field

#List all rows in intersection layer, for each row list polygon id, point value and object id
alldata = [row for row in arcpy.da.SearchCursor(in_table=fc, field_names=[polygon_id_field, point_value_field, 'OID@'])]
#[('66H 0h SV', 990.197, 653238), ('66H 0h SV', 983.55, 656320), ('66H 0h SV', 965.07, 651390), ('66H 0h SV', 964.652, 657665)] #Like this

alldata = sorted(alldata, key=lambda x: (x[0], x[1]), reverse=True) #Sort by polygon id and point value so the lowest value comes last

#('66H 0h SV', 990.197, 653238)
#('66H 0h SV', 983.55, 656320)

minval = {val[0]:val[2] for val in alldata} #For each polygon id, store the lowest points object id
#{'66H 0h SV': 650386, '66H 0g SV': 651881...     #So for polygon '66H 0h SV', point with object id 650386 has the lowest DOY value
to_select = tuple(minval.values())
#(650386, 651881, 657618, 653111, 650398, 655888, ...
sql = "{0} IN{1}".format(arcpy.AddFieldDelimiters(datasource=fc, field=arcpy.Describe(fc).OIDFieldName), to_select) #Create a sql clause
#'"OBJECTID" IN(650386, 651881, 657618, 653111, 650398, ...)

arcpy.SelectLayerByAttribute_management(in_layer_or_view=fc, where_clause=sql) #Select them

enter image description here

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.