Geometry : Forms of the same feature .... # 7

550
0
05-13-2019 08:25 AM
Labels (1)
DanPatterson_Retired
MVP Emeritus
2 0 550

Taking things apart and putting them back together isn't as easy as one would think

Consider the multipart shape with two holes.  The first part is the larger of the two.

The coordinates can be derived from the polygon in a variety of ways.  The table below shows it blown apart using the __geo_interface__ method.

Multipart polygon from __geo_interface__Converted to an array
s01
<Polygon object at 
   0x1e7faa2a128[0x1e7faa396c0]>

geo1 = s01.__geo_interface__['coordinates']

geo1
 
[[[(300010.0, 5000020.0),
   (300010.0, 5000010.0),
   (300000.0, 5000010.0),
   (300000.0, 5000020.0),
   (300010.0, 5000020.0)],
  [(300002.0, 5000018.0),
   (300002.0, 5000012.0),
   (300008.0, 5000012.0),
   (300008.0, 5000018.0),
   (300002.0, 5000018.0)]],
 [[(300007.0, 5000017.0),
   (300007.0, 5000013.0),
   (300003.0, 5000013.0),
   (300003.0, 5000017.0),
   (300007.0, 5000017.0)],
  [(300005.0, 5000016.0),
   (300004.0, 5000014.0),
   (300006.0, 5000014.0),
   (300005.0, 5000016.0)]]]‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍
a_min = [300000, 5000000]
geo1a = [(np.array(i) - a_min)
         for prt in geo1
         for i in prt]

geo1a
 
[array([[10., 20.],
        [10., 10.],
        [ 0., 10.],
        [ 0., 20.],
        [10., 20.]]),
 array([[ 2., 18.],
        [ 2., 12.],
        [ 8., 12.],
        [ 8., 18.],
        [ 2., 18.]]),
 array([[ 7., 17.],
        [ 7., 13.],
        [ 3., 13.],
        [ 3., 17.],
        [ 7., 17.]]),
 array([[ 5., 16.],
        [ 4., 14.],
        [ 6., 14.],
        [ 5., 16.]])]
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

We now have a basis of comparison.  Note that the polygon is converted to its multipart point form.

Points
shapes = arcpy.da.FeatureClassToNumPyArray(in_fc0,
               ['OID@', 'SHAPE@X', 'SHAPE@Y'],
               "", SR, True)

shapes
array([(2, 300010., 5000020.), (2, 300010., 5000010.),
       (2, 300000., 5000010.), (2, 300000., 5000020.),
       (2, 300010., 5000020.), (2, 300002., 5000018.),
       (2, 300002., 5000012.), (2, 300008., 5000012.),
       (2, 300008., 5000018.), (2, 300002., 5000018.),
       (2, 300007., 5000017.), (2, 300007., 5000013.),
       (2, 300003., 5000013.), (2, 300003., 5000017.),
       (2, 300007., 5000017.), (2, 300005., 5000016.),
       (2, 300004., 5000014.), (2, 300006., 5000014.),
       (2, 300005., 5000016.)],
      dtype=[('OID@', '<i4'),
             ('SHAPE@X', '<f8'),
             ('SHAPE@Y', '<f8')])

z = shapes[shapes['OID@'] == 2]
z0 = stu(z[['SHAPE@X', 'SHAPE@Y']]) - [300000, 5000000]‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Points was derived using FeatureClassToNumPyArray

The approach is outlined in the table to the right.

The final points list is given in the table below along with their representation in polygon, polyline and line form.

PolygonPolylineLinePoints
[Geo([[10., 20.],
      [10., 10.],
      [ 0., 10.],
      [ 0., 20.],
      [10., 20.],
      [nan, nan],
      [ 2., 18.],
      [ 2., 12.],
      [ 8., 12.],
      [ 8., 18.],
      [ 2., 18.]]),
 Geo([[ 7., 17.],
      [ 7., 13.],
      [ 3., 13.],
      [ 3., 17.],
      [ 7., 17.],
      [nan, nan],
      [ 5., 16.],
      [ 4., 14.],
      [ 6., 14.],
      [ 5., 16.]])]‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍
[Geo([[ 0., 10.],
      [10., 10.],
      [10., 20.],
      [ 0., 20.],
      [ 0., 10.]]),
 Geo([[ 5., 16.],
      [ 4., 14.],
      [ 6., 14.],
      [ 5., 16.]]),
 Geo([[ 7., 17.],
      [ 7., 13.],
      [ 3., 13.],
      [ 3., 17.],
      [ 7., 17.]]),
 Geo([[ 2., 18.],
      [ 2., 12.],
      [ 8., 12.],
      [ 8., 18.],
      [ 2., 18.]])]‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍
Geo([[ 0., 10.],
     [10., 10.],
     [10., 20.],
     [ 0., 20.],
     [ 0., 10.],
     [ 5., 16.],
     [ 4., 14.],
     [ 6., 14.],
     [ 5., 16.],
     [ 7., 17.],
     [ 7., 13.],
     [ 3., 13.],
     [ 3., 17.],
     [ 7., 17.],
     [ 2., 18.],
     [ 2., 12.],
     [ 8., 12.],
     [ 8., 18.],
     [ 2., 18.]])‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍
array([[10., 20.],
       [10., 10.],
       [ 0., 10.],
       [ 0., 20.],
       [10., 20.],
       [ 2., 18.],
       [ 2., 12.],
       [ 8., 12.],
       [ 8., 18.],
       [ 2., 18.],
       [ 7., 17.],
       [ 7., 13.],
       [ 3., 13.],
       [ 3., 17.],
       [ 7., 17.],
       [ 5., 16.],
       [ 4., 14.],
       [ 6., 14.],
       [ 5., 16.]])‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

--------------------------------------------

Putting them back together is hard.

If you have simple shapes, you can reassemble a point set using NumPyArrayToFeatureClass

If the shape is complex (multiple parts and/or holes) then you can't reassemble properly.

That will be the topic of the next blog

About the Author
Retired Geomatics Instructor at Carleton University. I am a forum MVP and Moderator. Current interests focus on python-based integration in GIS. See... Py... blog, my GeoNet blog...
Labels