3

I am interested in converting a polygon featureclass created via Create Fishnet (Data Management) from singleparts to multipart. QGIS has a very handy tool called singleparts to multipart that does the job, however, I need to stay within the arcpy/Python realm. The ArcGIS equivalent appears to be Dissolve, yet the tool does not preserve individual feature geometry when the polygons share a border.

The following is an example of the type of fishnet polygon data I am working with, which was generated from the below script.

enter image description here

import arcpy

# Create a fishnet with 9 columns and 9 rows
# with origin at (1, 1) and output geometry is set to default (POLYGON)
arcpy.CreateFishnet_management("C:/temp/fishnet2.shp", "1 1", "1 9", "1", "1", "9",     "9", "#", "NO_LABELS", geometry_type = "POLYGON")

How can I convert a fishnet polygon from singleparts to multipart using a scripting approach?

*Note. There is a similar question at GIS SE (Create a single connected polygon from multiple polygons), although it appears to be limited to desktop tools. I am looking for any Python or arcpy-based approaches.

Community
  • 1
Aaron
  • 51,658
  • 28
  • 154
  • 317
  • 3
    Also related: http://gis.stackexchange.com/questions/63042/how-to-create-touching-multi-part-polygons (although not specifically related to scripting). I'd be interested to hear if Menno's answer works! – Erica Oct 18 '14 at 23:01
  • I agree with @Erica - this sounds like a duplicate of the question she cited. – PolyGeo Oct 18 '14 at 23:09
  • @PolyGeo Please note that this is directed at any Python approach, which may or may not include the use of the arcpy site package. I have not found any questions that relate to a scripted Python approach. – Aaron Oct 18 '14 at 23:37
  • I just updated my answer at http://gis.stackexchange.com/a/63065/115 to clarify that I think this is a limitation of shapefiles and geodatabase feature classes i.e. it seems to be independent of whether they are written by a product/component from the ArcGIS platform, Python, or another GIS vendor or open source GIS project. – PolyGeo Oct 18 '14 at 23:57
  • @PolyGeo Seeing that you can convert adjacent (Fishnet polygons) singleparts to multipart in QGIS, this is likely not a limitation of the shapefile format--but rather, a limitation in out-of-the-box ArcGIS tools. I will test Menno's cursor approach and update the post... – Aaron Oct 19 '14 at 00:19
  • 1
    You're right - I had overlooked the significance of that sentence in your question and have not used that QGIS tool. I'm keen to learn the outcome of your test because, as I said in my answer to the other question, I could see no logic in the limitation being in place for the ArcGIS platform. If you are able to create the desired multipart polygon with internal boundaries using ArcPy I am curious to know whether it can be drawn in both QGIS and ArcGIS for Desktop, and whether it passes the Check Geometry tool test of ArcGIS for Desktop. – PolyGeo Oct 19 '14 at 00:31
  • 2
    While the QGIS tool may do it, this would seem to suggest the results are invalid geometries. I wonder, similar to PolyGeo, if the results pass a geometry check in QGIS, never mind in ArcGIS. I recall mention that one is more strict than the other when it comes to geometry. I think the question referenced in the question's Note is sufficiently different from what you want to do as that question is only creating a single, contiguous polygon. The one Erica references is much closer and was the first thing that came to my mind when I saw this question as well. – Chris W Oct 19 '14 at 02:00
  • I just tested the script myself, unfortunately it didn't work, see updated answer. – Menno Oct 20 '14 at 08:06
  • @PolyGeo QGIS is able to generate the multipart output, however, there are geometry errors as you and Chris W indicated. When the geometry is repaired, the vector features do not resemble the original data (see Menno's screenshots). – Aaron Oct 20 '14 at 12:41
  • Obviously based on this question and the linked 'touching multipart' question we can come up with reasons for doing this. I'm curious as to that rationale here - why would you want all the squares retained individually rather than just the outer boundary if what they represent is identical (ie single attribute row)? In the linked question it's clear there is a difference between one of the three parts and they just want to ignore it in terms of records/attributes but not geometry. – Chris W Oct 20 '14 at 19:11
  • @Chris W I was working on a tool to help visually estimate vegetative cover. The actual output would have been the intersection of fishnet polygons and forest stand polygons. The idea was to have clusters of multipart fishnet polygons that I could assign a total cluster count to the attribute table and then assign an estimated count to another field. The end result would have been for example 20 veg polygons/100 total polygons = 20% cover for cluster X. – Aaron Oct 20 '14 at 19:57
  • I see. I'd wondered about maybe using lines instead of polygons for the fishnet and making a multipart feature out of those if a fill or area wasn't critical to the use. – Chris W Oct 20 '14 at 21:22
  • That could be a good option @Chris W if it were feasible to count the number of polyline "cells" in each cluster. – Aaron Oct 20 '14 at 21:33
  • 1
    Create label points with the fishnet creation, get a point count, then discard the points as an intermediate layer? I confess I don't fully understand the process and products so this may be more problems than it solves. – Chris W Oct 20 '14 at 21:45

2 Answers2

4

I am not sure if the following will not do the same thing as dissolve, but if I'm correct, it should not.

You can use a SearchCursor() to loop through the polygons, get each polygon's geometry, add these as parts on a new polygon geometry object, and use an InsertCursor() to insert this new record. 

sc = arcpy.SearchCursor("c:/temp/fishnet2.shp")
ic = arcpy.InsertCursor("yournewfeatureclass.shp")
newRow = ic.newRow()
array = arcpy.Array()
for row in sc:
    geom = row.getValue("shape").getPart(0)
    array.add(geom)
polygon = arcpy.Polygon(array)
newRow.setValue("shape", polygon)
ic.insertRow(newRow)
del sc, ic

Note: this inserts all of the polygons if your initial feature class as one multipart polygon in the output feature class. You can use MakeFeatureLayer_management() and SelectLayerByAttribute_management() to select the right polygons to be added together.

Okay, so I didn't test the above and it didn't work, unfortunately (anyone know how to strikethrough text?). Below is the result of my test. Curiously, the effect wasn't the same for an FGDB feature class and a shapefile. The polygons in the searchcursor are sorted by Objectid, the bottomleft being the first, traveling first upward along colums, then along rows.

Effect of adding neighbouring polygons in 1 geometry

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
Menno
  • 590
  • 2
  • 11
  • 1
    Have you tested this? Shapefiles are more forgiving for topological errors, but they remain errors (in violation of the specification), and the data would be useless for analysis without repair. – Vince Oct 19 '14 at 16:21
  • I have now, see my updated answer. – Menno Oct 20 '14 at 08:01
  • 1
    ArcGIS tries to repair invalid polygon geometries; the difference is due to ring orientation differences. – Vince Oct 20 '14 at 10:50
  • Thanks Menno. These are the results I was getting after doing repair geometry on the QGIS output. – Aaron Oct 20 '14 at 12:37
  • 1
    @Aaron This is sort of what I expected. Consider a valid multipart that is a hole inside an outer boundary - it's clear what the hole is and what the boundary is. But in your desired case with multiple shared edges, how would it know what is a hole and what isn't? The desire is that none of them are holes, but an edge indicates something is changing. This reminds me of something said in another question (I think by Vince actually) that a polygon can't contain itself, which is exactly what you'd get with the hole situation if the shared edge didn't indicate change. – Chris W Oct 20 '14 at 19:04
3

I think that you will have to keep the polygons as separate features. Polygon multiparts are an evil construct that has no logical topology for analysis and even display in your case.

What is wrong with a separate featureclass? If you want to aggregate the properties, then a Dissolve will do this. You could relate the multiple features to the single dissolved outline if you want a single attribute record. Alternatively you could convert the polygons into polylines to remove the duplicate coincident rings. There are ways of building topology rules to relate the boundary lines and the polygon.

You need valid, unambiguous topology in any system to be able to do more analysis. Consider simple ideas such as a polygon centroid. It can fall outside all the polygons so a point-in polygon overlay will fail to select the set.

Coverages handle this case easily because polygons were defined using a PAL (Polygon-Arc-List) that defined a polygon as the chain of lines. No overlaps possible (unless you went to a more complex Region featureclass). I miss coverage topology, it made this sort of analysis a no-brainer.

kimo
  • 371
  • 2
  • 7