topic Re: Create a points grid for each polygon in a polygon file in Geoprocessing Questions

Create a points grid for each polygon in a polygon file

KevinCressy — Tue, 08 Jun 2010 10:33:12 GMT

Hello all,

I have found the following python script (Generate Lattice - http://arcscripts.esri.com/details.asp?dbid=14398) that will create a points grid.

I have run the script for a small area and it works very well. however I was wondering if anyone could advise on how I could modify the script for my purposes?

I have been given a polygon file that contains many polygons covering a very large area, however there is a much larger area that does not contain any polygons at all. I have been asked to add a point grid to each of the polygons. One way I thought of doing this would be to use the script and then select only the points that fall within the polygons - however the problem I have is that the extent of the file is so large I end up creating thousands of points that I don't need, plus this process of creating the initial grid takes many hours to run.

I can understand that the way to approach this would be as below, I'm just not sure about what steps I need to take to achieve this:

1) Go the first polygon in my polygon file
2) extract the min x and min y coords of that features envelope
3) create a points grid from the origin of the the min x and y coords
4) move to the next polygon.

Could anyone give any advice on how I might be able to achieve this?

Any help would be appreciated.

Regards,

Kevin

Re: Create a points grid for each polygon in a polygon file

CedricWannaz — Tue, 08 Jun 2010 16:49:43 GMT

Hello Kevin,

Are the extents of your polygons (not the polygons) overlapping or are they completely disjoint?

Regards,

Cedric

Re: Create a points grid for each polygon in a polygon file

CedricWannaz — Sat, 11 Dec 2021 18:21:48 GMT

If you are interested in an un-elegant solution, I spent 20 minutes and wrote what follows. It is certainly possible to do everything in a few smart operations using the Desktop interface, but I am not that smart, and elegance is not always mandatory 😉

The code produces the output that is attached to this post. I was not sure whether you wanted one single point feature class in the end or as many as there are polygons, so I didn't implement the union..

Cheers,

Cedric

def buildSubArray( fc, extent, dist_x, dist_y ):

    # Determine the number of points that can fit in the x and y directions
    # in the extent of the polygon, and create a centered grid in this extent.

    # Note: points are equidistant in the chosen coordinate system (and NOT
    # on ground, unless you are working with an equidistant - hence rather 
    # local - coordinate system).

    range_x   = extent.xmax-extent.xmin
    nPoints_x = (int) (range_x / dist_x)
    start_x   = extent.xmin + 0.5 * (range_x - nPoints_x * dist_x)

    range_y   = extent.ymax-extent.ymin
    nPoints_y = (int) (range_y / dist_y)
    start_y   = extent.ymin + 0.5 * (range_y - nPoints_y * dist_y)

    point  = gp.CreateObject( 'POINT' )
    cursor = gp.InsertCursor( fc )
    
    for i_y in range( nPoints_y ) :
        for i_x in range( nPoints_x ) :
            feature = cursor.NewRow()
            point.x = start_x + i_x * dist_x
            point.y = start_y + i_y * dist_y
            feature.shape = point
            cursor.InsertRow( feature )
    
    del feature
    del cursor
    del point
    


print '--- Start.\n'


# -- Initialization.

# Create and initialize the Geoprocessor.
import arcgisscripting

gp = arcgisscripting.create( 9.3 )
gp.SetProduct( 'ArcInfo' )
gp.Workspace = r'C:\Users\wannaz\Documents\Projects\ArcGIS General\CreateArrayOfPoints\project.gdb'
gp.OverwriteOutput = 1

# Define setup-specific parameters.
polygon_fc      = 'polygons'
temporary_fc    = '__points'
clipped_fc_base = 'clippedPoints_'

xSpacing        = 5
ySpacing        = xSpacing

# Create a layer from the polygon feature class for select() purpose.
polygon_layer = 'polygons_layer'
gp.MakeFeatureLayer_management( polygon_fc, polygon_layer )

# Extract the OBJECTID field name.
polygon_desc = gp.Describe( polygon_layer )
polygon_OIDFieldName = polygon_desc.OIDFieldName


# -- Loop over polygons.

rows = gp.SearchCursor( polygon_layer ) 
row  = rows.Next()

while row :
    oid = row.GetValue( polygon_OIDFieldName )    
    print ' Processing polygon', oid, '..',  

    # - Create array of points that covers the extent of current polygon.
    
    gp.CreateFeatureclass( gp.Workspace, temporary_fc, "POINT", polygon_fc, "", "", polygon_desc.SpatialReference )
    extent = row.shape.Extent
    buildSubArray( temporary_fc, extent, xSpacing, ySpacing )
    
    # - Select current polygon and clip.

    whereClause = polygon_OIDFieldName + " = " + str( oid )
    gp.SelectLayerByAttribute_management( polygon_layer, "NEW_SELECTION", whereClause ) 
    
    clipped_fc = clipped_fc_base + str( oid )
    gp.Clip_analysis( temporary_fc, polygon_layer, clipped_fc )

    print "done."
    row  = rows.Next()

del rows
del row    


# -- Union?

# If you were interested in the union of sub grids but not in the sub grids
# themselves, you could create them in_memory to be faster. In such case, you
# could create a list of sub grid names in the loop above, that you would pass
# to the Union_analysis() tool/method. You would have to delete the temporary
# sub grids afterwards.  
# PS: it would take about two minutes for me to update this script for it to
#    perform the union, so let me know if it would help.


# -- Mr. Proper.

gp.Delete_management( temporary_fc )
gp.Delete_management( polygon_layer )


print '\n--- End of Job.'

Re: Create a points grid for each polygon in a polygon file

KevinCressy — Thu, 10 Jun 2010 12:31:51 GMT

Hi Cedric,

Unfortunately it appears that your solution requires an ArcInfo licence - I only have ArcView so have been unable to test this.

Re: Create a points grid for each polygon in a polygon file

CedricWannaz — Thu, 10 Jun 2010 19:21:04 GMT

Did you try to run it after replacing "ArcInfo" with "ArcView" in the SetProduct() call ?

Re: Create a points grid for each polygon in a polygon file

KevinCressy — Mon, 14 Jun 2010 10:17:32 GMT

Hello all, not sure if someone can help here, I seem to be getting quite close now - but I am still struggling.

I have put together the following script to create a grid of points within the boundaries of each polygon in my polygon file.

I am however getting an unexpected error in Python - which is telling me that it cannot find the created shapefile even though I see that these are created in the output directory. The error I am getting is below:

ERROR 000210: Cannot create output C:\Python_Test\outputBELppend_fishnet_pnt.shp
Failed to execute (Clip).

What seems strange is that the file path is not correct in the returned error. The file path should be:

C:\Python_Test\output\append_fishnet_pnt.shp

and NOT

C:\Python_Test\outputBELppend_fishnet_pnt.shp

i.e. the work BEL seems to have replaced '\a'

Any idea why this might be?

The full script is given below:

import arcgisscripting

gp = arcgisscripting.create(9.3)
gp.overwriteOutput = 1

inpolys = r'C:\Python_Test\Polygons.shp'

infc = inpolys
#dsc = gp.Describe(infc)
#sr = dsc.spatialreference
#gp.outputCoordinateSystem = sr

rows = gp.SearchCursor(infc)

try:
row = rows.next() # take me to first row

id = 0 # this index is used to rename the output for each iteration

while row:
feat = row.shape
ext = feat.Extent
xmin = ext.XMin
xmax = ext.XMax
ymin = ext.YMin
ymax = ext.YMax
origin = str(xmin) + " " + str(ymin)
direction = ymin + float(100)
opp_corner = str(xmax) + " " + str(ymax)
y_axis = str(xmin) + " " + str(direction)
outfc = r"C:\Python_Test\output\myfishnettemp.shp" # + str(i)
gp.CreateFishnet_management(outfc, origin, y_axis, "0.5", "0.5", "0", "0", opp_corner, 'LABELS')

#fid = 0 # as it is a shapefile, it starts with zero for row in rows: # you are using just to iterate not to use the row object
outIndivPolygonfc = r"C:\Python_Test\output\indiv_poly.shp"
where_clause = '"FID" = %i' % (id)
gp.Select_analysis(infc, outIndivPolygonfc, where_clause)

gp.Clip_analysis("C:\Python_Test\output\myfishnettemp_label.shp", "C:\Python_Test\output\indiv_poly.shp", "C:\Python_Test\output\append_fishnet_pnt.shp")

gp.Append_management('C:\Python_Test\output\append_fishnet_pnt.shp', "C:\Python_Test\output\final_fishnet_pnt.shp", "NO_TEST")

row = rows.next()

id = id + 1

except:
print gp.GetMessages(2)

Re: Create a points grid for each polygon in a polygon file

CedricWannaz — Mon, 14 Jun 2010 12:22:48 GMT

Correct:
input = r"C:\aProject\newFile.shp"
input = "C:\\aProject\\newFile.shp"
input = "C:/aProject/newFile.shp"

Incorrect:
input = "C:\aProject\newFile.shp"

Using single or double quotes doesn't change the behavior.

Note: I wrote what follows initially but this is wrong; I didn't see your Select():
"""By the way, if the extent of a feature/polygon overlaps another polygon in your setup, your method will fail because the clip will keep overlapping points. This is why I asked about the extents on 06-08-2010 09:49 AM, and why the example code that I wrote above "isolates" each polygon before clipping""".

Cedric

Re: Create a points grid for each polygon in a polygon file

KevinCressy — Tue, 15 Jun 2010 05:16:13 GMT

I eventually managed to crack this with help from Peter McDaid and Nobbir Ahmed, Craig and Cedric. Thank you everyone else who replied to my threads and emails!

I include my final code below. Note that I have also added additional functionality to the script which now creates square buffers from the point features.

import arcgisscripting

gp = arcgisscripting.create(9.3)
gp.overwriteOutput = 1

# Example Argument Inputs: C:\Python_Test\output\final_fishnet_pnt.shp C:\Python_Test\output\final_fishnet_sqr.shp 0.5 TRUE
inPoints = sys.argv[1] #Input point feature class - must be shapefile
outPolys = sys.argv[2] #Output polygon feature class - must be shapefile
bufDist = float(sys.argv[3]) #Buffer distance
keepFields = sys.argv[4] #Boolean type: Maintain fields and field values of the input in the output

#This is the input polygon
inpolys = r"C:\Python_Test\Polygons.shp"

infc = inpolys

rows = gp.SearchCursor(infc)

try:
# take me to first row in polygon file
row = rows.next()

# this index is used as the fid value in the where clause query
id = 0

while row:
feat = row.shape
ext = feat.Extent
xmin = ext.XMin
xmax = ext.XMax
ymin = ext.YMin
ymax = ext.YMax
origin = str(xmin) + " " + str(ymin)
direction = ymin + float(100)
opp_corner = str(xmax) + " " + str(ymax)
y_axis = str(xmin) + " " + str(direction)
outfc = r"C:\Python_Test\output\myfishnettemp.shp"
gp.CreateFishnet_management(outfc, origin, y_axis, "5", "5", "0", "0", opp_corner, 'LABELS')

#fid = 0 # as it is a shapefile, it starts with zero for row in rows: # you are using just to iterate not to use the row object
outIndivPolygonfc = r"C:\Python_Test\output\indiv_poly.shp"
where_clause = '"FID" = %i' % (id)
gp.Select_analysis(infc, outIndivPolygonfc, where_clause)

gp.Clip_analysis(r"C:\Python_Test\output\myfishnettemp_label.shp", r"C:\Python_Test\output\indiv_poly.shp", r"C:\Python_Test\output\append_fishnet_pnt.shp")

gp.Append_management(r"C:\Python_Test\output\append_fishnet_pnt.shp", r"C:\Python_Test\output\final_fishnet_pnt.shp", "NO_TEST")

row = rows.next()

id = id + 1

# Prepare the output based on whether field and field values are desired in the output
if keepFields:
# Create empty output polygon feature class that includes fields of the input
gp.CreateFeatureClass(os.path.dirname(outPolys), os.path.basename(outPolys), "POLYGON", inPoints, "", "", inPoints)

# Create a short-list of fields to ignore when moving fields values from input to output
ignoreFields = []
# Use Describe properties to identify the ShapeFieldName and OIDFieldName
desc = gp.Describe(inPoints)
ignoreFields.append(desc.ShapeFieldName)
ignoreFields.append(desc.OIDFieldName)

# Create a list of fields to use when moving field values from input to output
fields = gp.ListFields(inPoints)
#field = fields.Next()
#fieldList = []
#while field:
# if field.Name not in ignoreFields:
# fieldList.append(field.Name)
# field = fields.Next()

fieldList = []
for field in fields:
if field.Name not in ignoreFields:
fieldList.append(field.Name)

else:
# Create empty output polygon feature class without fields of the input
gp.CreateFeatureClass(os.path.dirname(outPolys), os.path.basename(outPolys), "POLYGON", "", "", "", inPoints)

# Open searchcursor
inRows = gp.SearchCursor(inPoints)
inRow = inRows.Next()

# Open insertcursor
outRows = gp.InsertCursor(outPolys)

# Create point and array objects
pntObj = gp.CreateObject("Point")
arrayObj = gp.CreateObject("Array")

while inRow: # One output feature for each input point feature
inShape = inRow.Shape
pnt = inShape.GetPart(0)
# Need 5 vertices for square buffer: upper-right, upper-left, lower-left,
# lower-right, upper-right. Add and subtract distance from coordinates of
# input point as appropriate.
for vertex in [0,1,2,3,4]:
pntObj.id = vertex
if vertex in [0,3,4]:
pntObj.x = pnt.x + bufDist
else:
pntObj.x = pnt.x - bufDist
if vertex in [0,1,5]:
pntObj.y = pnt.y + bufDist
else:
pntObj.y = pnt.y - bufDist
arrayObj.add(pntObj)

# Create new row for output feature
feat = outRows.NewRow()

# Shift atttributes from input to output
if keepFields:
for fieldName in fieldList:
feat.SetValue(fieldName, inRow.GetValue(fieldName))

# Assign array of points to output feature
feat.Shape = arrayObj
# Insert the feature
outRows.InsertRow(feat)

# Clear array of points
arrayObj.RemoveAll()

# Get next feature in searchcursor
inRow = inRows.Next()

# Delete inputcursor
del outRows

except:
print gp.GetMessages(2)

Re: Create a points grid for each polygon in a polygon file

IanLane — Mon, 31 Jul 2023 20:03:51 GMT

This is a pretty old post, but if you're still around I wouldn't mind seeing the union code. I was able to implement your code with some tweaks