Select to view content in your preferred language

Dispersing geocoded points within a Polygon

6542
15
Jump to solution
11-27-2013 03:56 AM
MTP
by
New Contributor II
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)])
Tags (2)
0 Kudos
15 Replies
JohnShell
New Contributor II

Xander

Thank you for your reply.

I still get an error message for now.

Thank you.

John

Error below:

Traceback (most recent call last):

  File "H:\DispersePoints.py", line 49, in <module>

    main()

  File "H:\DispersePoints.py", line 22, in main

    disperse_points(fcPoints,polygon)

  File "H:\DispersePoints.py", line 40, in disperse_points

    while not inside:

UnboundLocalError: local variable 'inside' referenced before assignment

 

My code:

 

#------------------------------------------------------------------------------- 

# Name:        Disperse3.py 

# Purpose:     Disperse points in multiple polygons 

# Author:      arcpy Team  xander_bakker  https://community.esri.com/people/xander_bakker

#              http://arcpy.wordpress.com/2013/06/07/disperse-overlapping-points/ 

# Created:     02-dec-2013  

#------------------------------------------------------------------------------- 

 

import arcpy, random

 

def main(): 

       

    fcPoints = r"H:\\CondoAddressPoints.gdb\\new_points2"  

    fcPolygon = r"H:\\CondoAddressPoints.gdb\\new_polys" 

    

  

    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]):

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

0 Kudos
DarrenWiens2
MVP Honored Contributor

Try initializing 'inside' to None just inside the p loop:

for p in points:
     inside = None‍‍‍‍

edit: compare your indentation to Xander's at 'while not inside'. It's different and is responsible for your problem.

JohnShell
New Contributor II

After making the change in indentation the script runs perfectly.

Thank you everyone

John

0 Kudos
MTP
by
New Contributor II
Hello Xander,

thank you very much. I will try to get it running today. As soon as I have the results I post them here. It is exactly what I was searching for.

Best regards
Martin
0 Kudos
MTP
by
New Contributor II
Hello Xander,

the script is working perfectly. I even managed to run it with 500.000 points.

Here is an example of a smaller region.

[ATTACH=CONFIG]29614[/ATTACH]

Thank you very much
Martin
0 Kudos
XanderBakker
Esri Esteemed Contributor
Hi Martin,

I'm glad I could be of any assistance. If you think it was helpful you    can use the "arrow" button in order to help other  members find useful   information:



More info here: http://resources.arcgis.com/en/help/forums-mvp/

Kind regards,

Xander
0 Kudos