[arcpy] how to merge polygons into one multipart polygon

309
3
Jump to solution
09-25-2020 05:09 PM
VincentLantaca
New Contributor II

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!

Tags (1)
0 Kudos
1 Solution

Accepted Solutions
DanPatterson
MVP Honored Contributor
 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.


... sort of retired...

View solution in original post

3 Replies
DanPatterson
MVP Honored Contributor
 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.


... sort of retired...

View solution in original post

VincentLantaca
New Contributor II

Thank you!

0 Kudos
JoshuaBixby
MVP Esteemed Contributor

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"
)