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; } } }

Create a points grid for each polygon in a polygon file

KevinCressy

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

Vind meer berichten getagd met

Accepted answers

All comments

CedricWannaz

Hello Kevin,

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

Regards,

Cedric

CedricWannaz

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.'

Screenshot.pdf

project.gdb.rar

KevinCressy

Hi Cedric,

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

CedricWannaz

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

KevinCressy

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)

CedricWannaz

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

KevinCressy

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)

IanLane

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

Doe mee aan het gesprek

Aangemelde leden kunnen berichten plaatsen, updates volgen en meer. Nieuw hier? Registreer een gratis account.

Snelle Links

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