Clip FishNet Grids by Geometry Shape

1241
8
Jump to solution
01-30-2019 03:39 PM
AdamThompson3
New Contributor III

Hey Everyone,

Currently I am using the tool found here;

http://www.arcgis.com/home/item.html?id=9398bd2232cb4c8490b0b05015364d28  

Now it works great and does exactly what id like but one slight tweak I would like to make is as it iterates each shape and creates the grid, I would just like to Intersect it with my original grids geometry so I do not get any grids outside the shape. For example; 

Now I believe I have identified the area of the code in which I need to perform this operation but it has not been working for me, wondering if anyone could point me in the direction I am possibly going wrong or help in anyway possible. Here is the area of code I am changing and inserting my tweaks. 

The error returned is shown below; 

Thank you to anyone who has the time to help out! Appreciate it. 

0 Kudos
1 Solution

Accepted Solutions
RandyBurton
MVP Alum

I've done some experimenting with the code you found for creating polygon fishnets.  It adds some steps to clip the fishnet to the shape of the original polygons.  I did not explore the group options in the code, so there may need to be some adjustments for these options to work correctly in the modified code.

The basic changes involve tracking the object IDs of the starting polygons and clipping the fishnet boxes using the geometry intersect method.

In my test, I had 3 polygons that I wanted to cover with a 20x20 grid and clip that to the shapes of the polygons.  The results looked like:

polygons with clipped fishnet

Here's the code starting around line 63 (following the set-up and functions in the code linked to in your question):

# this would be line 63 following the functions and other set-up code
# ## read in_features geometry into dictionary for use with geometry intersect method later
origPolys = { r[0] : r[1] for r in arcpy.da.SearchCursor(in_features,['OID@','SHAPE@']) } # ## list comprehension

# Process the data
scratchName = arcpy.CreateScratchName("MBG","ByWidth","FeatureClass")
arcpy.MinimumBoundingGeometry_management(in_features,scratchName,"RECTANGLE_BY_WIDTH",group_option,group_fields,"MBG_FIELDS")
arcpy.DeleteField_management(scratchName,["MBG_Width","MBG_Length"]) # ## ["MBG_Width","MBG_Length","ORIG_FID"]) remove ORIG_ID from delete field list
arcpy.CreateFeatureclass_management(os.path.dirname(out_features),\
                                    os.path.basename(out_features),"POLYGON",\
                                    template=scratchName,\
                                    spatial_reference=sR)
cur = arcpy.da.InsertCursor(out_features,["SHAPE@","MBG_Orientation","ORIG_FID"]) # ## ["SHAPE@","MBG_Orientation"]) add ORIG_ID to field list
if labels == "LABELS":
    arcpy.CreateFeatureclass_management(os.path.dirname(out_features),\
                                        os.path.basename(out_features)+"_label","POINT",\
                                        spatial_reference=sR)
for row in arcpy.da.SearchCursor(scratchName,["OID@","SHAPE@","MBG_Orientation","ORIG_FID"]): # ## ["OID@","SHAPE@","MBG_Orientation"]): add ORIG_ID to field list
    oid = row[0]
    ofid = row[3] # ## original feature id ( in_features objectid )
    arcpy.AddMessage("\nProcessing polygon " + str(oid))
    mbgPoly = row[1]
    orientation = float(row[2])
    pivotPoint = mbgPoly.firstPoint
    origin_coord = str(pivotPoint.X) + ' ' + str(pivotPoint.Y)
    orthogonalPoly = RotatePolygon(mbgPoly,pivotPoint,orientation)
    opposite_corner = str(orthogonalPoly.extent.XMax) + ' ' + str(orthogonalPoly.extent.YMax)
    yCoord = str(pivotPoint.X) + ' ' + str(pivotPoint.Y + 1.0)
    scratchName2 = arcpy.CreateScratchName("OrthogonalFishNet",str(oid),"FeatureClass")
    arcpy.CreateFishnet_management(scratchName2,origin_coord,yCoord,cell_width,cell_height,number_rows,number_columns,labels=labels,template=orthogonalPoly.extent,geometry_type="POLYGON")
    #arcpy.CreateFishnet_management(scratchName2,origin_coord,yCoord,cell_width,cell_height,number_rows,number_columns,opposite_corner,labels=labels,geometry_type="POLYGON")

    # ## create a buffer from the original polygon which is in the dictionary
    polyBuff = origPolys[ofid].buffer(0) # ## buffer distance either 0 or negative number in default spatial reference units
    # arcpy.AddMessage("\npolyBuff = {}".format(polyBuff.JSON)) # ## for debugging
   
    for cell in arcpy.da.SearchCursor(scratchName2,["OID@","SHAPE@"]):
        arcpy.AddMessage("Processing fishnet cell..."  + str(cell[0]))
        outPoly = RotatePolygon(cell[1],pivotPoint,0-orientation)
        
        # arcpy.AddMessage("\noutPoly = {}".format(outPoly.JSON)) # ## for debugging (json shows spatial reference as null)
        # ## since outPoly has a null spatial reference, we need to specify it for the intersect method
        newPoly = arcpy.FromWKT(outPoly.WKT,sR) # ## use WKT and add default spatial reference http://desktop.arcgis.com/en/arcmap/latest/analyze/arcpy-functions/fromwkt.htm
        # arcpy.AddMessage("\nnewPoly = {}".format(newPoly.JSON)) # ## for debugging
        # ## clip outPoly as required using geometry's intersect method
        bufPoly = newPoly.intersect(polyBuff,4) # ## http://desktop.arcgis.com/en/arcmap/latest/analyze/arcpy-classes/geometry.htm  (see intersect method)
        # arcpy.AddMessage("\nbufPoly = {}".format(bufPoly.WKT)) # ## for debugging
        # ## it is possible that the polygon falls outside the intersect clip
        # ## json will show: bufPoly = {"rings":[],"spatialReference":{"wkid":102100,"latestWkid":3857}} and WKT will indicate EMPTY
        if "EMPTY" not in bufPoly.WKT: # ## add only if geometry is not empty (null)
            cur.insertRow([bufPoly,orientation,ofid])# ## [outPoly,orientation]) added field to save original polygon's objectid, see insertcursor around line 75

    arcpy.Delete_management(scratchName2)
    if labels == "LABELS":
        del cur #Otherwise two InsertCursors in the same workspace requires arcpy.da.Editor
        lcur = arcpy.da.InsertCursor(out_features+"_label",["SHAPE@"])
        for label in arcpy.da.SearchCursor(scratchName2+"_label",["SHAPE@"]):
            lcur.insertRow([arcpy.PointGeometry(RotatePoint(label[0].firstPoint,pivotPoint,0-orientation),sR)])
        arcpy.Delete_management(scratchName2+"_label")
        del lcur
        cur = arcpy.da.InsertCursor(out_features,["SHAPE@","MBG_Orientation"])
arcpy.Delete_management(scratchName)
del cur
arcpy.AddMessage("\nDone")‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Hope this helps.

View solution in original post

8 Replies
DanPatterson_Retired
MVP Emeritus

polygons are expected for both the inputs, is your red line a polygon or polyline?

0 Kudos
AdamThompson3
New Contributor III

It is a polygon just with a hollow fill for visual purposes, but I actually do believe I see why the error is occuring. Within the code scratchName refers to the Minimum Bounding Box created, and it is iterating through each bounding box created instead of the original polygons shape. So essentially I am trying to intersect the rotated fishnet with the bounding box it created, and I believe there lies the problem.

I am trying to now contain this all under a search cursor in which it iterates through each original input polygon shape instead of each bounding box shape, that way I can call on the feature it is iterating through. Example of what I am playing around with; 

0 Kudos
RandyBurton
MVP Alum

I am assuming you wish to remove the part of the fishnet that is outside the red lines.

The section of code you highlighted is a possibility for adding the code for clipping.  The issue I see is that a link to the original polygon is lost in the process by the DeleteField_management "ORIG_FID" line in the code.  There are several places that the ORIG_FID would need to be referenced in some of the lines that follow.  The reference to mbgPoly is to the minimum bounding box polygon, so it wouldn't do any clipping to the fishnet.

The use of multiple cursors in that section of code would also be problematic as noted in the nearby comment.  It might be possible to read the required geometry into a dictionary that could be referenced.  And, I'm sure there are other solutions, too.

0 Kudos
AdamThompson3
New Contributor III

Yea that is what I am trying to achieve

I will look into the dictionary solution, that is a good suggestion. Currently outside of that script I have been able to achieve what Id like with this code; 

just a matter of finding out how to merge that within the code. 

0 Kudos
RandyBurton
MVP Alum

I've done some experimenting with the code you found for creating polygon fishnets.  It adds some steps to clip the fishnet to the shape of the original polygons.  I did not explore the group options in the code, so there may need to be some adjustments for these options to work correctly in the modified code.

The basic changes involve tracking the object IDs of the starting polygons and clipping the fishnet boxes using the geometry intersect method.

In my test, I had 3 polygons that I wanted to cover with a 20x20 grid and clip that to the shapes of the polygons.  The results looked like:

polygons with clipped fishnet

Here's the code starting around line 63 (following the set-up and functions in the code linked to in your question):

# this would be line 63 following the functions and other set-up code
# ## read in_features geometry into dictionary for use with geometry intersect method later
origPolys = { r[0] : r[1] for r in arcpy.da.SearchCursor(in_features,['OID@','SHAPE@']) } # ## list comprehension

# Process the data
scratchName = arcpy.CreateScratchName("MBG","ByWidth","FeatureClass")
arcpy.MinimumBoundingGeometry_management(in_features,scratchName,"RECTANGLE_BY_WIDTH",group_option,group_fields,"MBG_FIELDS")
arcpy.DeleteField_management(scratchName,["MBG_Width","MBG_Length"]) # ## ["MBG_Width","MBG_Length","ORIG_FID"]) remove ORIG_ID from delete field list
arcpy.CreateFeatureclass_management(os.path.dirname(out_features),\
                                    os.path.basename(out_features),"POLYGON",\
                                    template=scratchName,\
                                    spatial_reference=sR)
cur = arcpy.da.InsertCursor(out_features,["SHAPE@","MBG_Orientation","ORIG_FID"]) # ## ["SHAPE@","MBG_Orientation"]) add ORIG_ID to field list
if labels == "LABELS":
    arcpy.CreateFeatureclass_management(os.path.dirname(out_features),\
                                        os.path.basename(out_features)+"_label","POINT",\
                                        spatial_reference=sR)
for row in arcpy.da.SearchCursor(scratchName,["OID@","SHAPE@","MBG_Orientation","ORIG_FID"]): # ## ["OID@","SHAPE@","MBG_Orientation"]): add ORIG_ID to field list
    oid = row[0]
    ofid = row[3] # ## original feature id ( in_features objectid )
    arcpy.AddMessage("\nProcessing polygon " + str(oid))
    mbgPoly = row[1]
    orientation = float(row[2])
    pivotPoint = mbgPoly.firstPoint
    origin_coord = str(pivotPoint.X) + ' ' + str(pivotPoint.Y)
    orthogonalPoly = RotatePolygon(mbgPoly,pivotPoint,orientation)
    opposite_corner = str(orthogonalPoly.extent.XMax) + ' ' + str(orthogonalPoly.extent.YMax)
    yCoord = str(pivotPoint.X) + ' ' + str(pivotPoint.Y + 1.0)
    scratchName2 = arcpy.CreateScratchName("OrthogonalFishNet",str(oid),"FeatureClass")
    arcpy.CreateFishnet_management(scratchName2,origin_coord,yCoord,cell_width,cell_height,number_rows,number_columns,labels=labels,template=orthogonalPoly.extent,geometry_type="POLYGON")
    #arcpy.CreateFishnet_management(scratchName2,origin_coord,yCoord,cell_width,cell_height,number_rows,number_columns,opposite_corner,labels=labels,geometry_type="POLYGON")

    # ## create a buffer from the original polygon which is in the dictionary
    polyBuff = origPolys[ofid].buffer(0) # ## buffer distance either 0 or negative number in default spatial reference units
    # arcpy.AddMessage("\npolyBuff = {}".format(polyBuff.JSON)) # ## for debugging
   
    for cell in arcpy.da.SearchCursor(scratchName2,["OID@","SHAPE@"]):
        arcpy.AddMessage("Processing fishnet cell..."  + str(cell[0]))
        outPoly = RotatePolygon(cell[1],pivotPoint,0-orientation)
        
        # arcpy.AddMessage("\noutPoly = {}".format(outPoly.JSON)) # ## for debugging (json shows spatial reference as null)
        # ## since outPoly has a null spatial reference, we need to specify it for the intersect method
        newPoly = arcpy.FromWKT(outPoly.WKT,sR) # ## use WKT and add default spatial reference http://desktop.arcgis.com/en/arcmap/latest/analyze/arcpy-functions/fromwkt.htm
        # arcpy.AddMessage("\nnewPoly = {}".format(newPoly.JSON)) # ## for debugging
        # ## clip outPoly as required using geometry's intersect method
        bufPoly = newPoly.intersect(polyBuff,4) # ## http://desktop.arcgis.com/en/arcmap/latest/analyze/arcpy-classes/geometry.htm  (see intersect method)
        # arcpy.AddMessage("\nbufPoly = {}".format(bufPoly.WKT)) # ## for debugging
        # ## it is possible that the polygon falls outside the intersect clip
        # ## json will show: bufPoly = {"rings":[],"spatialReference":{"wkid":102100,"latestWkid":3857}} and WKT will indicate EMPTY
        if "EMPTY" not in bufPoly.WKT: # ## add only if geometry is not empty (null)
            cur.insertRow([bufPoly,orientation,ofid])# ## [outPoly,orientation]) added field to save original polygon's objectid, see insertcursor around line 75

    arcpy.Delete_management(scratchName2)
    if labels == "LABELS":
        del cur #Otherwise two InsertCursors in the same workspace requires arcpy.da.Editor
        lcur = arcpy.da.InsertCursor(out_features+"_label",["SHAPE@"])
        for label in arcpy.da.SearchCursor(scratchName2+"_label",["SHAPE@"]):
            lcur.insertRow([arcpy.PointGeometry(RotatePoint(label[0].firstPoint,pivotPoint,0-orientation),sR)])
        arcpy.Delete_management(scratchName2+"_label")
        del lcur
        cur = arcpy.da.InsertCursor(out_features,["SHAPE@","MBG_Orientation"])
arcpy.Delete_management(scratchName)
del cur
arcpy.AddMessage("\nDone")‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Hope this helps.

AdamThompson3
New Contributor III

Amazing! Thanks Randy for looking into that!

0 Kudos
BruceHarold
Esri Regular Contributor
AdamThompson3
New Contributor III

Thanks Bruce! I honestly had no idea this tool existed  but this is awesome as well! Thank you for the link 

0 Kudos