How to create centroid of points within a certain radius?

BinCai · ‎11-16-2010

We have ~30k points and points close to each other (for example 50 meters or less) should be grouped together as one point (centroid), as shown in attached screenshot (green dots are points we have and blue dots are centroids we want). ET Geowizards can do it but we the demo version has limit of 100 features. Is there any other free tool or method can do it? Thanks!

Bin

ChrisSnyder · ‎01-20-2011

Here it is.... Like I indicated earlier, this is a pretty crude way to do this, and I'm sure there are better ways now in v10.

# Description
# -----------
# This script assigns a unique identifier to features that are connected/close to each other.
# This may be useful for identifying interconnected networks of polygons or lines such as a road
# or stream networks. It is particularly useful when trying to do large single part dissolves...
# just dissolve by the unique identifier field (or in addition to other fields as well).
#
# Written By: Chris Snyder, WA DNR, 7/27/2007, chris.snyder(at)wadnr.gov
#
# Written For: Python 2.4.1 and ArcGIS v9.2 SP2
#
#Notes on input parameters:
#VARIABLE           PAREMETER_INDEX     PARAMETER_DATA_TYPE
#-------------------------------------------------------------------
#inLayer                 0                   Feature Layer
#networkIdField          1                   Text (for example, "NET_ID") #new long int field that will be added
#selectionMethod         2                   Text (for example, "INTERSECT") - this needs to be one of the selection keywords used in the SelectByLocation tool
#distThreshold           3                   Linear Unit (for example, "1 METER") only used if selectionMethod = "WITHIN_A_DISTANCE"

# Import system modules and creates the Geoprocessor object
try:
    
    #Import system modules and creates the Geoprocessor object (the v9.2 way)
    import sys, string, os, glob, shutil, time, traceback, arcgisscripting
    gp = arcgisscripting.create()
        
    #Defines some functions used for getting messages from the gp and python
    def showGpMessage():
        gp.AddMessage(gp.GetMessages())
        print >> open(logFile, 'a'), gp.GetMessages()
        print gp.GetMessages()
    def showGpWarning():
        gp.AddWarning(gp.GetMessages())
        print >> open(logFile, 'a'), gp.GetMessages()
        print gp.GetMessages()
    def showGpError():
        gp.AddError(gp.GetMessages())
        print >> open(logFile, 'a'), gp.GetMessages()
        print gp.GetMessages()
    def showPyMessage():
        gp.AddMessage(str(time.ctime()) + " - " + message)
        print >> open(logFile, 'a'), str(time.ctime()) + " - " + message
        print str(time.ctime()) + " - " + message
    def showPyWarning():
        gp.AddWarning(str(time.ctime()) + " - " + message)
        print >> open(logFile, 'a'), str(time.ctime()) + " - " + message
        print str(time.ctime()) + " - " + message
    def showPyError():
        gp.AddError(str(time.ctime()) + " - " + message)
        print >> open(logFile, 'a'), str(time.ctime()) + " - " + message
        print str(time.ctime()) + " - " + message

    #Sets a date/time stamp variable in the format yyyymmddhhmmss (example: 20060711132345)
    dateTimeStamp = time.strftime('%Y%m%d%H%M%S')

    #Specifies the root directory variable, defines the logFile variable, and does some minor error checking...
    dateTimeString = str(time.strftime('%Y%m%d%H%M%S'))
    scriptName = os.path.split(sys.argv[0])[-1].split(".")[0]
    userName = string.lower(os.environ.get("USERNAME")).replace(" ","_").replace(".","_")
    tempPathDir = os.environ["TEMP"]
    logFileDirectory = r"\\snarf\am\div_lm\ds\gis\tools\log_files"
    if os.path.exists(logFileDirectory) == True:
        logFile = os.path.join(logFileDirectory, scriptName + "_" + userName + "_" + dateTimeString + ".txt")
        try:
            print >> open(logFile, 'a'), "Write test successfull!"
        except:
            logFile = os.path.join(tempPathDir, scriptName + "_" + userName + "_" + dateTimeString + ".txt")  
    else:
        logFile = os.path.join(tempPathDir, scriptName + "_" + userName + "_" + dateTimeString + ".txt")
    if os.path.exists(logFile)== True:
        os.remove(logFile)
        message = "Created log file " + logFile; showPyMessage()
    message = "Running " + sys.argv[0]; showPyMessage()
    
    #Attempts to check out the lowest grade license available...
    if gp.CheckProduct("ArcView") == "Available":
        gp.SetProduct("ArcView")
    elif gp.CheckProduct("ArcEditor") == "Available":
        gp.SetProduct("ArcEditor")
    elif gp.CheckProduct("ArcInfo") == "Available":
        gp.SetProduct("ArcInfo")
    message =  "Selected an " + gp.ProductInfo() + " license"; showPyMessage()
    
    message = "Starting geoprocessing routine..."; showPyMessage()
    #*****************GEOPROCESSING STUFF GOES HERE******************************

    #Process: Sets some gp settings
    gp.overwriteoutput = 1

    #Defines the input feature layer  
    inLayer = gp.GetParameterAsText(0)
    networkIdField = gp.GetParameterAsText(1)
    selectionMethod = gp.GetParameterAsText(2)
    distThreshold = gp.GetParameterAsText(3)

    #Process: Tests for the existance of networkIdField - if it doesn't exists, add it
    if gp.listfields(inLayer, networkIdField).next() == None:
        gp.AddField_management(inLayer, networkIdField, "LONG"); showGpMessage()
    else:
        message = "Specified existing field type is: " + str(gp.listfields(inLayer, networkIdField).next().type); showPyMessage() 
        if gp.listfields(inLayer, networkIdField).next().type not in ["SmallInteger","Double","Integer"]:
            message = "ERROR: Existing specified network id field type must be Long Integer, Short Integer or Double! Exiting script..."; showPyError(); sys.exit()
        gp.CalculateField_management(inLayer, networkIdField, "0", "VB"); showGpMessage() #makes sure the field is nulled

    #Process: Figures out how many features are in the input feature class
    featureCount = int(gp.getcount(gp.describe(inLayer).catalogpath))

    #Process: Figures out what the oidValue field is called
    oidFieldName = gp.describe(inLayer).oidfieldname

    #Process: Does a recursive select by location
    featureLayer1 = "feature_layer_1"    
    featureLayer2 = "feature_layer_2"
    identifiedFeatures = 0
    networkId = 0
    gp.MakeFeatureLayer_management(inLayer, featureLayer1, "")
    while identifiedFeatures < featureCount:
        
        #Process: Assign a counter
        networkId =  networkId + 1
        
        #Process: Establishes a cursor in order to return a single OID of a non-identified feature
        try: #this is in a try loop cause' it will fail on the last loop since networkIdField will be (hopefully) 100% populated
            oidValue = gp.searchcursor(inLayer, networkIdField + " = 0 OR " + networkIdField + " IS NULL").next().getvalue(oidFieldName)
        except:
            message = "Breaking selection loop..."; showPyMessage()
            break #break out of the loop
        gp.SelectLayerByAttribute_management(featureLayer1, "NEW_SELECTION", oidFieldName + " = " + str(oidValue))

        #Process: Keep selecting features until there are no more to select
        count1 = 0
        count2 = 1
        count3 = 0
        while count2 > count1:
            count3 = count3 + 1 
            count1 = int(gp.getcount(featureLayer1))
            if selectionMethod == "WITHIN_A_DISTANCE":
                gp.SelectLayerByLocation_management(featureLayer1, "WITHIN_A_DISTANCE", featureLayer1, distThreshold, "NEW_SELECTION")
            else:
                gp.SelectLayerByLocation_management(featureLayer1, selectionMethod, featureLayer1, "", "NEW_SELECTION")
            count2 = int(gp.getcount(featureLayer1))
            if divmod(count3, 10)[1] == 0:
                message = "Processed the " + str(count3) + "th loop of connected feature ID #" + str(networkId) + "..."; showPyMessage()
        gp.CalculateField_management(featureLayer1, networkIdField, networkId, "VB")

        #Process: Every 25th loop, spit out a status message (too many status messages makes the script run slow) 
        if divmod(networkId, 25)[1] == 0:
            gp.MakeFeatureLayer_management(inLayer, featureLayer2, networkIdField + " > 0")
            identifiedFeatures = int(gp.getcount(featureLayer2))
            message = "Assigned identifier to " + str(float(identifiedFeatures) / float(featureCount) * 100) + " percent of the features"; showPyMessage()    

    #Process: Gives a final report
    gp.MakeFeatureLayer_management(inLayer, featureLayer2, networkIdField + " > 0")
    identifiedFeatures = int(gp.getcount(featureLayer2))
    gp.MakeFeatureLayer_management(inLayer, featureLayer2, networkIdField + " = 0 OR " + networkIdField + " IS NULL")
    nonIdentifiedFeatures = int(gp.getcount(featureLayer2))
    message = "CONNECTED FEATURES REPORT:"; showPyWarning()
    message = "--------------------------"; showPyWarning()
    message = "Assigned " + str(networkId -1) + " unique network identifiers to " + str(identifiedFeatures) + " features!"; showPyWarning()
        
    #*****************GEOPROCESSING STUFF ENDS HERE******************************
    message = "Ending geoprocessing routine..."; showPyMessage()

    #Indicates that the script is complete - both a visual and physical indicator...
    message = sys.argv[0] + " is all done!"; showPyMessage()
except:
    message = "\n*** LAST GEOPROCESSOR MESSAGE (may not be source of the error)***"; showPyError()
    showGpMessage()
    message = "PYTHON TRACEBACK INFO: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyError()
    message = "PYTHON ERROR INFO: " +  str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyError()