[arcpy] how to merge polygons into one multipart polygon

[VincentLantaca]

I am trying to merge polygons into one multi part polygon. I am trying to just append the polygon parts into one array and then creating a polygon from that array:

import arcpy

target_layer = r'D:\DatabaseConnection\SCRAPWORK.sde\egdb.TEST.TEST_V'

where = "OBJECTID IN (1601,1602,1603)"

arcpy.MakeFeatureLayer_management(target_layer, 'test_layer', where_clause=where)

polygons = []

arr = []
sr = None

with arcpy.da.SearchCursor('test_layer', ['SHAPE@']) as cursor:
    for row in cursor:
        polygon = row[0]

        print('isMultipart: {}'.format(polygon.isMultipart))
        print('partCount: {}'.format(polygon.partCount))
        sr = polygon.spatialReference
        print('')

        # add the polygon to an array
        arr.append(polygon.getPart())

array = arcpy.Array(arr)
multiPartPolygon = arcpy.Polygon(array, sr)

print('MULTI PART POLYGON:')
print('JSON:\n{}'.format(multiPartPolygon.JSON))
print('isMultipart: {}'.format(multiPartPolygon.isMultipart))
print('partCount: {}'.format(multiPartPolygon.partCount))
print('spatialReference: {}'.format(multiPartPolygon.spatialReference))
print('getPart: {}'.format(multiPartPolygon.getPart()))‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

When I run the script, I am seeing that the three features that I am trying to merge are only 1 part polygons, but the polygon created has 0 part count, and no data in rings:

isMultipart: False
partCount: 1

isMultipart: False
partCount: 1

isMultipart: False
partCount: 1

MULTI PART POLYGON:
JSON:
{"rings":[],"spatialReference":{"wkid":4326,"latestWkid":4326}}
isMultipart: True
partCount: 0
spatialReference: <geoprocessing spatial reference object object at 0x1065FC80>
getPart: <geoprocessing array object object at 0x1065FC80>‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Can anyone point out what I am doing wrong? Any help is greatly appreciated!

[DanPatterson · Accepted solution]

 arr = []
polys = []
sr = "NAD 1983 CSRS MTM  9"
with arcpy.da.SearchCursor('test', ['SHAPE@']) as cursor:
    for row in cursor:
        polygon = row[0]
        for cnt in range(polygon.partCount):
            part = polygon.getPart(cnt)
            polys.append(part)
            arr.append(part)
a = arcpy.Array(arr)

a
<Array 
  [<Array [<Point (300009.0, 5000001.0, #, #)>, <Point (300000.0, 5000000.0, #, #)>,
           <Point (300002.0, 5000008.0, #, #)>, <Point (300008.0, 5000010.0, #, #)>,
           <Point (300010.0, 5000010.0, #, #)>, <Point (300010.0, 5000008.0, #, #)>,
           <Point (300009.0, 5000001.0, #, #)>,
           None,
          <Point (300003.0, 5000003.0, #, #)>, <Point (300007.0, 5000003.0, #, #)>,
          <Point (300005.0, 5000007.0, #, #)>, <Point (300003.0, 5000003.0, #, #)>]>,
   <Array [<Point (300010.0, 5000008.0, #, #)>, <Point (300010.0, 5000010.0, #, #)>,
           <Point (300008.0, 5000010.0, #, #)>, <Point (300008.0, 5000011.0, #, #)>,
           <Point (300008.0, 5000012.0, #, #)>, <Point (300012.0, 5000012.0, #, #)>,
           <Point (300012.0, 5000008.0, #, #)>, <Point (300010.0, 5000008.0, #, #)>]>,
   <Array [<Point (300008.0, 5000011.0, #, #)>, <Point (300005.0, 5000010.0, #, #)>,
           <Point (300005.0, 5000012.0, #, #)>, <Point (300006.0, 5000012.0, #, #)>,
           <Point (300008.0, 5000012.0, #, #)>, <Point (300008.0, 5000011.0, #, #)>]>,
   <Array [<Point (300006.0, 5000012.0, #, #)>, <Point (300005.0, 5000012.0, #, #)>,
           <Point (300005.0, 5000015.0, #, #)>, <Point (300007.0, 5000014.0, #, #)>,
           <Point (300006.0, 5000012.0, #, #)>]>]>
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Before and after.  polygon with hole and 3 separate polygons.

The output (green)

p0 = arcpy.Polygon(a)

p0.partCount

p0.isMultipart
True‍‍‍‍‍‍

p0.getPart(0)
<Array [<Point (300006.0001220703, 5000012.00012207, #, #)>,
        <Point (300008.0001220703, 5000012.00012207, #, #)>,
        <Point (300012.0001220703, 5000012.00012207,
        <Point (300006.0001220703, 5000012.00012207, #, #)>,
        None,
        <Point (300003.0001220703, 5000003.00012207, #, #)>,
        <Point (300007.0001220703, 5000003.00012207, #, #)>,
        <Point (300005.0001220703, 5000007.00012207, #, #)>,
        <Point (300003.0001220703, 5000003.00012207, #, #)>]>‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Unexpected Result!!!!

Not at all, since you can't create multiparts that are disjoint if they share an edge, so the polygon resultant is indeed a multipart shape since it contains a hole, but the perimeter of the green output is the result of the dissolved boundaries from the input singlepart polygons.

[VincentLantaca]

Thank you!

[JoshuaBixby]

I have found that unioning creates more expected results, especially if there is a chance of polygons overlapping.  For example, looking at the following three single-part polygons:

import arcpy
from arcpy.management import CopyFeatures

SR = arcpy.SpatialReference(26915)
polys_wkt = [
    "POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))",
    "POLYGON((5 5, 15 5, 15 15, 5 15, 5 5))",
    "POLYGON((0 20, 10 20, 10 30, 0 30, 0 20))"
]

CopyFeatures(
    [arcpy.FromWKT(wkt, SR) for wkt in polys_wkt],
    "memory/polys"
)
‍‍‍‍‍‍‍‍‍‍‍‍‍‍

What do you expect the new multi-part polygon to look like?  This

CopyFeatures(
    arcpy.Polygon(
        arcpy.Array([arcpy.FromWKT(wkt, SR).getPart(0) for wkt in polys_wkt]),
        SR
    ),
    "memory/polys_array"
)
‍‍‍‍‍‍‍

or this?

from functools import reduce

CopyFeatures(
    reduce(arcpy.Geometry.union,[arcpy.FromWKT(wkt, SR) for wkt in polys_wkt]),
    "memory/polys_union"
)
‍‍‍‍‍‍