I'm currently working with a point shapefile derived from OSM that displays intersecting vertices of a road network. This shapefile has had duplicate geometries removed, resulting in a single point per link intersection. Although OSM is rather close to actual road topology, this leaves a junction with divided roads using 4 points at 1 road junction. I need to merge these 4 points into 1 to create a single point at a junction. There are occurrences where the junction will have 3 points, these will also need to be merged to 1 center point as well. A graphic displaying what needs to be done is attached. The blue points are the original junction points, while the green center point is what is needed. The shapefile consists of over 20,000 unique points, therefore manually completing this task is not efficient.


I have used Arcmap 10.6 and ArcGIS Pro (similar tools) to try converting the roads into polygons and getting the center line from the polygon, where I can then split at intersections and get single points at intersections, however this didn't work for extracting the correct points. I have also tried snapping points within a distance, however this didn't create the new center point I'm looking for. I have also tried dissolving the point shapefile based on common attributes, but the attributes are not consistent enough and this also doesn't create a new center point. I have spent hours researching this topic, but always came up empty handed.

ArcMap 10.6 and ArcGIS Pro with advanced license. My scripting knowledge is almost non-existent.

  • Does the four Points in your example have some group attribute that is the same four all four Points?
    – BERA
    Aug 22, 2018 at 6:37
  • There are many that experience similar attributes, such as an identical highway tag, however this is not consistent throughout the file and may experience two different road types intersecting. The only merging criteria I see as common throughout the file is 1 point's spatial relationship to the other 3. Aug 22, 2018 at 6:41
  • Please decide which of QGIS, ArcGIS Pro and ArcMap you wish to ask about in this particular question. That way you can describe precisely what you have tried and where you are stuck with that.
    – PolyGeo
    Aug 22, 2018 at 15:49

1 Answer 1


You can use arcpy. Center Point will only be created when number of Points in a Group is between 2 and 4. Execute code in the python window of Desktop/Pro after changing the four indicated lines:

import arcpy
from itertools import combinations

#Adjust the following four lines:
arcpy.env.workspace = r'C:\Test' #Can also be a geodatabase, if so remove .shp below
fc = r'Junctions.shp' #Input fc located in workspace
new_fc = r'Junction_centers.shp' #Output fc, will be created in workspace
groupingdistance = 50 #Group points together that are within this distance of eachother.

#Group points by proximity
def make_equiv_classes(pairs): #https://stackoverflow.com/questions/39915402/combine-a-list-of-pairs-tuples
    groups = {}
    for (x, y) in pairs:
        xset = groups.get(x, set([x]))
        yset = groups.get(y, set([y]))
        jset = xset | yset
        for z in jset:
            groups[z] = jset
    return set(map(tuple, groups.values()))

all_points = [i for i in arcpy.da.SearchCursor(fc,['OID@','SHAPE@'])]
groups = []
for combo in combinations(all_points,2):
    if combo[0][1].distanceTo(combo[1][1]) <= groupingdistance:
groups = make_equiv_classes(sorted(groups, key=lambda x: (x[0],x[1])))

#Calculate mean x and y for each group and create a center point
arcpy.CreateFeatureclass_management(out_path=arcpy.env.workspace, out_name=new_fc, 
icur = arcpy.da.InsertCursor(new_fc,'SHAPE@XY')
for junction in groups:
    if 1<len(junction)<5: #Only create center points if number of points in a group is between 2 and 4
        sql = """{0} IN({1})""".format(arcpy.AddFieldDelimiters(fc,arcpy.Describe(fc).OIDFieldName),','.join([str(i) for i in junction]))
        xys = [i[0] for i in arcpy.da.SearchCursor(fc,'SHAPE@XY',sql)]
        meanx = sum([x[0] for x in xys])/len([x[0] for x in xys])
        meany = sum([y[1] for y in xys])/len([y[1] for y in xys])
del icur

enter image description here


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.