nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

How to create centroid of points within a certain radius?

BinCai

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

centroid.jpg

Amanda__Huber_1-1747676768816.png

Find more posts tagged with

Accepted answers

All comments

Hornbydd

Bin,

If ET Geowizard works for you, why don't you buy it and support it's continued development?

Also have a look at the Measuring Geographic Distributions toolset in the Spatial Statistics toolbox.

Duncan

BinCai

Thanks, Duncan. My boss doesn't want to pay for it. 😞
I will check the measuring toolset.

DanPatterson_Retired

How accurately do you need the centroid?
Are the clusters identified with a particular attribute? If so, then a bounding container of some form (ie extent rectangle, min area bounding rectangle, convex hull), followed by a centroid or label point determination might be reasonable.
If the clusters have no identified attribute which can separate the clusters, then perhaps a dissolved buffer of some reasonable distance to ensure polygons with no holes, followed by a similar centrality point determination might suffice.
otherwise, you may have to resort to other clustering algorithms (k-means?) etc.
In any event, your 30k points may have to be batched into separate geographic regions followed by a merge.
If your example is not really representative of the actual patterning in your data, then all bets are off until further information can be provided

LaurenRosenshein

Hi Bin,

Another option for this type of grouping is to start with the Integrate tool (Data Management), which allows you to snap features together that are within a specified distance of each other. The integrate tool will change the coordinates of all input features so that those features within a specified distance of each other will be coincident. Please note that the Integrate tool will change your input feature class, so it is always a good idea to run Copy Features beforehand to ensure that you do not lose your original data. Once you have run Integrate, you can use the Collect Events tool from the Spatial Statistics toolbox to create a count of the number of coincident points at each location.

This workflow is detailed in the Hot Spot Analysis Tutorial as a way to aggregate incident data for a hot spot analysis. The tutorial provides step-by-step instructions on creating a model to Integrate and then Collect Events using point data.

Hope this helps!

Lauren Rosenshein
Geoprocessing Product Engineer

Check out the latest Spatial Statistics resources at http://bit.ly/spatialstats!

LaurenRosenshein

That is a great question, Bill! So great, in fact, that I had to go ask a colleague down the hall. 🙂 The response that I got is as follows:

A given point can move up to 1.4 X the cluster tolerance, but it�??s an iterative process, so the movement of a point, could move it within the cluster tolerance of other points, and shift it further. I think in this case, yes, all the points will shift, and will shift to the center of the grouping (2.5, 0).

My colleague also pointed me to the help topic, Topology in ArcGIS, which I think you will find provides a lot of additional information on cluster tolerances.

Hope this helps!

Lauren

ChrisSnyder

I wrote a Python script that would get you part way there: http://forums.arcgis.com/threads/17809-polygon-in-polygon-analysis-count-if?p=56528&viewfull=1#post56528

The script basically does a recursive SelectByLocation (maybe specify 1 mile selction radius?). The FC then gets a CLUSTER_ID field (you can name it whatever you want) added and populated. So if there were three "clusters" of points that were > 1 mile appart (but the points in the cluster were all < 1 mile appart) you would end up with 3 different CLUSTER_ID values. Then all you need to do is add the X and Y coordinates to the table, run the Summary Statistics tool (using CLUSTER_ID as the case field, and getting the mean X and Y coordinate values). Build a new FC based on teh output of the summary stats table. Presto.

DanPatterson_Retired

Given Bill's interesting example, what is returned...and what should be considered a cluster given such a circumstance?

AdrianRingenbach

I wrote a Python script that would get you part way there: http://forums.arcgis.com/threads/17809-polygon-in-polygon-analysis-count-if?p=56528&viewfull=1#post56528

The script basically does a recursive SelectByLocation (maybe specify 1 mile selction radius?). The FC then gets a CLUSTER_ID field (you can name it whatever you want) added and populated. So if there were three "clusters" of points that were > 1 mile appart (but the points in the cluster were all < 1 mile appart) you would end up with 3 different CLUSTER_ID values. Then all you need to do is add the X and Y coordinates to the table, run the Summary Statistics tool (using CLUSTER_ID as the case field, and getting the mean X and Y coordinate values). Build a new FC based on teh output of the summary stats table. Presto.

Hi csny490

I've got a similar problem as it is mentioned here in this Thread and I think, your python script would solve it perfectly! unfortunately, it's not available anymore under the link you pasted here. Could you upload it again? Or could you sent it to me personally?

Thanks a lot!

rinadr

AdrianRingenbach

csny490;57040 wrote:
The script basically does a recursive SelectByLocation [..]. The FC then gets a CLUSTER_ID field (you can name it whatever you want) added and populated.

To be more precise: I dont nead a centroid Point as above mentioned, I need a function or a script, that populates a new field in the Attribute table with a ID which is different for every cluster.

Thanks a lot!

ChrisSnyder

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()

Join the Conversation

Signed in members can post, follow updates, and more. New here? Register a free account.

Quick Links

Community Help
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.