Group nearby points and find center point for each group

Question

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.

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.append([combo[0][0],combo[1][0]])
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, 
                                    geometry_type='POINT', 
                                    spatial_reference=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])
        icur.insertRow([(meanx,meany)])
del icur