Dispersing geocoded points within a Polygon

[MTP]

I have posted the following question within the geocoding forum:
http://forums.arcgis.com/threads/97761-Dispersing-geocoded-points-within-a-Polygon

But perhaps it is more suitable to be listed here as it is more a Python question. Please use the link above for answers to prevent double posting. Please excuse the circumstances.

Original post:

Hi,

I have a problem with dispersing geocoded points (randomely/)equally within a polygon. More specific I have geocoded adresses using zip codes and would like to disperse them within a polygon. I know that it is possible to disperse markers but I would like to disperse the points. Dispersing them randomely would be not a problem as it is i.e. possible with Geowizard.

Therefore, I have searched the forums and have found some solutions for Avenue but in my ArcGIS 10 this obviously is not working. Based on the post here: http://arcpy.wordpress.com/2013/06/07/disperse-overlapping-points/#comment-220 I have tried to develop a script. But it is not working.

To visualize what I would like to achieve:

The geocoding result (with only one ZIP/normally I have around 500):
[ATTACH=CONFIG]29401[/ATTACH]

And the result I would like to receive:

[ATTACH=CONFIG]29402[/ATTACH]

Can anyone advice what could be changed to make it work? Is it necessary to have the x, y already in the table defined as columns?

Thank you very much!

import random import arcpy  arcpy.env.workspace = ???D:\???\database.gdb???  in_points = ???D:\???\???\geodatabase.gdb\Geocoding_Results.shp??? polygon = ???D:\???\geodatabase.gdb\zipcode.shp???  def point_in_poly(poly, x, y): ??????"Returns if the point is inside the polygon.  Parameters: poly: arcpy.Polygon() geometry x: x coordinate (float) y: y coordinate (float)  ??????" pg = arcpy.PointGeometry(arcpy.Point(x, y), poly.spatialReference) return poly.contains(pg)  def disperse_points(in_points, polygon): ??????"Randomly disperse points inside a polygon.  Parameters: in_points: Point feature class/layer (with or without selection) polygon: arcpy.Polygon() geometry  ??????"  lenx = polygon.extent.width leny = polygon.extent.height  with arcpy.da.UpdateCursor(in_points, "shape@xy") as points: for p in points: x = (random.random() * lenx) + polygon.extent.XMin y = (random.random() * leny) + polygon.extent.YMin inside = point_in_poly(polygon, x, y) while not inside: x = (random.random() * lenx) + polygon.extent.XMin y = (random.random() * leny) + polygon.extent.YMin inside = point_in_poly(polygon, x, y) points.updateRow([(x, y)])

[XanderBakker]

I tested with this code and it works as I would expect. The point is that you need to provide a polygon geometry and not a polygon featureclass to the disperse_points function.

#-------------------------------------------------------------------------------
# Name:        Disperse.py
# Purpose:     Disperse points in polygon
# Author:      arcpy Team
#              http://arcpy.wordpress.com/2013/06/07/disperse-overlapping-points/
# Created:     27-11-2013
#-------------------------------------------------------------------------------

def main():
    global arcpy, random
    import arcpy, random
    fcPoints = r"C:\Project\_Forums\Disperse\test.gdb\points" # edit this
    fcPolygon = r"C:\Project\_Forums\Disperse\test.gdb\polygon" # edit this

    with arcpy.da.SearchCursor(fcPolygon, "SHAPE@") as cursor:
        for row in cursor:
            polygon = row[0] # take first polygon
            break
    del row

    disperse_points(fcPoints,polygon)

    print "ready..."

def point_in_poly(poly, x, y):
    """Returns if the point is inside the polygon.

    Parameters:
        poly: arcpy.Polygon() geometry
        x:    x coordinate (float)
        y:    y coordinate (float)

    """
    pg = arcpy.PointGeometry(arcpy.Point(x, y), poly.spatialReference)
    return poly.contains(pg)

def disperse_points(in_points, polygon):
    """Randomly disperse points inside a polygon.

    Parameters:
        in_points:  Point feature class/layer (with or without selection)
        polygon:    arcpy.Polygon() geometry

    """

    lenx = polygon.extent.width
    leny = polygon.extent.height

    with arcpy.da.UpdateCursor(in_points, "SHAPE@XY") as points:
        for p in points:
            x = (random.random() * lenx) + polygon.extent.XMin
            y = (random.random() * leny) + polygon.extent.YMin
            inside = point_in_poly(polygon, x, y)
            while not inside:
                x = (random.random() * lenx) + polygon.extent.XMin
                y = (random.random() * leny) + polygon.extent.YMin
                inside = point_in_poly(polygon, x, y)
            points.updateRow([(x, y)])


if __name__ == '__main__':
    main()

Kind regards,

Xander

[MTP]

Thank you very much! Now the script is running. However all points get unique coordinates outside any of the zip polygons. It looks like this:

These are the overlapping points in the Zip areas:
[ATTACH=CONFIG]29415[/ATTACH]

After running the script I get this:
[ATTACH=CONFIG]29416[/ATTACH]

Both shapefiles have the same coordinate system.

Perhaps my error is here:

polygon = row[0] # take first polygon

Should this value be changed? In the ZIP file there are only the Polygons you see in the first picture.

In any case I have to thank you very much.

Best regards
MP

[RichardFairhurst]

Thank you very much! Now the script is running. However all points get unique coordinates outside any of the zip polygons. It looks like this:

These are the overlapping points in the Zip areas:
[ATTACH=CONFIG]29415[/ATTACH]

After running the script I get this:
[ATTACH=CONFIG]29416[/ATTACH]

Both shapefiles have the same coordinate system.

Perhaps my error is here:

polygon = row[0] # take first polygon

Should this value be changed? In the ZIP file there are only the Polygons you see in the first picture.

In any case I have to thank you very much.

Best regards
MP

I tested the script with a single polygon and 68 points and it ran perfectly.  It even worked when I used a multi-part polygon.  My projection was State Plane for both.  However, I only used 1 polygon and one set of points that fell inside of it.

Processing multiple polygons won't work if you use the raw feature class the way this is written.  You need to use layers and run the script for each polygon separately using a selection of the one polygon with its corresponding points.  Or rewrite the script to process all of the polygons with a test before moving a point that verifies the initial position of the point is inside the polygon and only dispersing those points that are inside.

[XanderBakker]

The points getting placed far away from any ZIP code polygon may be due to differences in coordinate system. Make sure your points and ZIP code polygons share the exact same coordinate system. The point coordinates are calculated using the extent of a single polygon.

Using multiple polygons would require some changes. You would have to loop through all the polygons and create a feature layer based on a select by location in which only those points are passed through that are within the current polygon. Another option would be a spatial join from the points to the polygons. By this you will probably obtain the FIPS at the points. Use this FIPS in a whereclause to obtain the right polygon and also in the whereclause to process only those points within that polygon.

Kind regards,

Xander

[MTP]

Thank you very much for your help. If I have a solution available I will post it here.

[XanderBakker]

To lend a hand:

• Perform a spatial join from your points to your polygons


• I assume your polygon fc has a numeric field which is unique for each polygon


• Define this field as the fldID

This code will loop over each polygon, read out the geometry and set up a where clause expression. This is fed to the disperse function where only the points in the polygon are dispersed.

This is a result of the code:
[ATTACH=CONFIG]29470[/ATTACH]


Another option is to change the disperse function and initially determine if the punt is within the current polygon and process the point only if it is inside (don't affect the point location if it's not within the current polygon)...

#-------------------------------------------------------------------------------
# Name:        Disperse2.py
# Purpose:     Disperse points in multiple polygons
# Author:      arcpy Team
#              http://arcpy.wordpress.com/2013/06/07/disperse-overlapping-points/
# Created:     29-11-2013
#-------------------------------------------------------------------------------

def main():
    global arcpy, random
    import arcpy, random
    fcPoints = r"C:\Project\_Forums\Disperse\test.gdb\points_joined" # points after spatial join
    fcPolygon = r"C:\Project\_Forums\Disperse\test.gdb\Polygons"
    fldID = 'myID' # polygon ID field in point FC after spatial join

    arcpy.env.overwriteOutput = True

    with arcpy.da.SearchCursor(fcPolygon, ("SHAPE@",fldID)) as cursor:
        for row in cursor:
            polygon = row[0]
            myId = row[1]
            expression = arcpy.AddFieldDelimiters(fcPoints, fldID) + " = {0}".format(myId)
            disperse_points(fcPoints,polygon,expression)
    del row
    print "ready..."

def point_in_poly(poly, x, y):
    pg = arcpy.PointGeometry(arcpy.Point(x, y), poly.spatialReference)
    return poly.contains(pg)

def disperse_points(in_points, polygon, expression):
    lenx = polygon.extent.width
    leny = polygon.extent.height

    with arcpy.da.UpdateCursor(in_points, "SHAPE@XY", where_clause=expression) as points:
        for p in points:
            x = (random.random() * lenx) + polygon.extent.XMin
            y = (random.random() * leny) + polygon.extent.YMin
            inside = point_in_poly(polygon, x, y)
            while not inside:
                x = (random.random() * lenx) + polygon.extent.XMin
                y = (random.random() * leny) + polygon.extent.YMin
                inside = point_in_poly(polygon, x, y)
            points.updateRow([(x, y)])


if __name__ == '__main__':
    main()

Kind regards,

Xander

[XanderBakker · Accepted solution]

What I mentioned in my previous post is that it isn't necessary to peform a spatial join on the points and polygons first. The code becomes more simple too.

By verifying first of a point is originally within the current polygon, you can process only those points within the current polygon.
This is done by the line:

            if point_in_poly(polygon, p[0][0], p[0][1]): 

The final code would be:

#------------------------------------------------------------------------------- # Name:        Disperse3.py # Purpose:     Disperse points in multiple polygons # Author:      arcpy Team #              http://arcpy.wordpress.com/2013/06/07/disperse-overlapping-points/ # Created:     02-dec-2013 #-------------------------------------------------------------------------------  def main():     global arcpy, random     import arcpy, random     fcPoints = r"C:\Project\_Forums\Disperse\test.gdb\points3"     fcPolygon = r"C:\Project\_Forums\Disperse\test.gdb\Polygons"      arcpy.env.overwriteOutput = True      with arcpy.da.SearchCursor(fcPolygon, ("SHAPE@")) as cursor:         for row in cursor:             polygon = row[0]             disperse_points(fcPoints,polygon)     del row     print "ready..."  def point_in_poly(poly, x, y):     pg = arcpy.PointGeometry(arcpy.Point(x, y), poly.spatialReference)     return poly.contains(pg)  def disperse_points(in_points, polygon):     lenx = polygon.extent.width     leny = polygon.extent.height      with arcpy.da.UpdateCursor(in_points, "SHAPE@XY") as points:         for p in points:             if point_in_poly(polygon, p[0][0], p[0][1]): # I changed code here!                 x = (random.random() * lenx) + polygon.extent.XMin                 y = (random.random() * leny) + polygon.extent.YMin                 inside = point_in_poly(polygon, x, y)                 while not inside:                     x = (random.random() * lenx) + polygon.extent.XMin                     y = (random.random() * leny) + polygon.extent.YMin                     inside = point_in_poly(polygon, x, y)                 points.updateRow([(x, y)])             else:                 pass # don't update location if point doesn't originally falls inside current polygon  if __name__ == '__main__':     main()

Kind regards,

Xander

[JohnShell]

Xander:

  I am fairly new to python and arcpy.

I was curious how to get your code Disperse3.py into the correct syntax and indentation where it will run correctly.

Thank you.

John Shell 

[XanderBakker]

Hi John, 

The code would probably be something like this:

#-------------------------------------------------------------------------------
# Name:        Disperse3.py
# Purpose:     Disperse points in multiple polygons
# Author:      arcpy Team
#              http://arcpy.wordpress.com/2013/06/07/disperse-overlapping-points/
# Created:     02-dec-2013
#-------------------------------------------------------------------------------

import arcpy
import random

def main():
    fcPoints = r"C:\Project\_Forums\Disperse\test.gdb\points3"
    fcPolygon = r"C:\Project\_Forums\Disperse\test.gdb\Polygons"
    arcpy.env.overwriteOutput = True

    with arcpy.da.SearchCursor(fcPolygon, ("SHAPE@")) as cursor:
        for row in cursor:
            polygon = row[0]
            disperse_points(fcPoints, polygon)
        del row
    print "ready..."

def point_in_poly(poly, x, y):
    pg = arcpy.PointGeometry(arcpy.Point(x, y), poly.spatialReference)
    return poly.contains(pg)

def disperse_points(in_points, polygon):
    lenx = polygon.extent.width
    leny = polygon.extent.height
    with arcpy.da.UpdateCursor(in_points, "SHAPE@XY") as points:
        for p in points:
            if point_in_poly(polygon, p[0][0], p[0][1]):
                # I changed code here!
                x = (random.random() * lenx) + polygon.extent.XMin
                y = (random.random() * leny) + polygon.extent.YMin
                inside = point_in_poly(polygon, x, y)
                while not inside:
                    x = (random.random() * lenx) + polygon.extent.XMin
                    y = (random.random() * leny) + polygon.extent.YMin
                    inside = point_in_poly(polygon, x, y)
                points.updateRow([(x, y)])
            else:
                pass # don't update location if point doesn't originally falls inside current polygon

if __name__ == '__main__':
    main()‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Just make sure to input points and polygons featureclasses mentioned on line 13 and 14 to point to the correct locations. Also, it would be better to work on a copy of the original data, since the code will update the point featureclass.