Skip navigation
All People > Dan_Patterson > Py... blog > 2019 > May
2019
Dan_Patterson

Geometry : Part 8

Posted by Dan_Patterson Champion May 30, 2019

Geometry

 

 

Geometry in NumPy... # 1

Geometry ... ArcPy and NumPy... # 2

Geometry ... Deconstructing poly* features # 3

Geometry ... Reconstructing Poly Features # 4

Geometry ... Attributes actually... the other bits # 5

Geometry: Don't believe what you see ... # 6

Geometry : Forms of the same feature # 7

 

Two multipart shapes, (one with holes, one without) and a singlepart shape.

 

Shapes

 

The shape itself is represented by an arcpy array of points.  Since there are two parts to the shape, there are two arrays.  The inner rings/holes are separated by None.

 

Each point carries the extra baggage of Z and M values whether they are needed or not.

 

 

 

Shape 1...

Object array... an array of arrays

Shapes 2 and 3...

ndarray...ndim = 3 shape = (2, 9, 2)

ndarray...ndim = 2 shape =  (4, 2)

Geo array (last 2 olumns, X, Y)

first 3 columns are for printing

array([
array([[10., 20.],  # first
       [10., 10.],  # outer
       [ 0., 10.],
       [ 0., 20.],
       [10., 20.],
       [nan, nan],
       [ 3., 19.],  # first
       [ 3., 13.],  # inner
       [ 9., 13.],
       [ 9., 19.],
       [ 3., 19.]]),
array([[ 8., 18.],  # second
       [ 8., 14.],  # outer
       [ 4., 14.],
       [ 4., 18.],
       [ 8., 18.],
       [nan, nan],
       [ 6., 17.],  # second
       [ 5., 15.],  # inner
       [ 7., 15.],
       [ 6., 17.]])],
dtype=object)

 

 

[<Array 
[<Point (300010.0, 5000020.0, #, #)>,
<Point (300010.0, 5000010.0, #, #)>,
<Point (300000.0, 5000010.0, #, #)>,
<Point (300000.0, 5000020.0, #, #)>,
<Point (300010.0, 5000020.0, #, #)>,
None,
<Point (300003.0, 5000019.0, #, #)>,
<Point (300003.0, 5000013.0, #, #)>,
<Point (300009.0, 5000013.0, #, #)>,
<Point (300009.0, 5000019.0, #, #)>,
<Point (300003.0, 5000019.0, #, #)>]>,
<Array
[<Point (300008.0, 5000018.0, #, #)>,
<Point (300008.0, 5000014.0, #, #)>,
<Point (300004.0, 5000014.0, #, #)>,
<Point (300004.0, 5000018.0, #, #)>,
<Point (300008.0, 5000018.0, #, #)>,
None,
<Point (300006.0, 5000017.0, #, #)>,
<Point (300005.0, 5000015.0, #, #)>,
<Point (300007.0, 5000015.0, #, #)>,
<Point (300006.0, 5000017.0, #, #)>]>]
array([[[12., 18.], # first
        [12., 12.],
        [20., 12.],
        [20., 10.],
        [10., 10.],
        [10., 20.],
        [20., 20.],
        [20., 18.],
        [12., 18.]],

       [[25., 24.],  # second
        [25., 14.],
        [15., 14.],
        [15., 16.],
        [23., 16.],
        [23., 22.],
        [15., 22.],
        [15., 24.],
        [25., 24.]]],
      dtype('float64'))

 

 

 

array([[14., 20.],
       [10., 20.],
       [15., 28.],
       [14., 20.]])
pnt shape  part  X       Y     
--------------------------------
000     0         10.00   20.00
001     0         10.00   10.00
002     0          0.00   10.00
003     0          0.00   20.00
004     0         10.00   20.00
005     0   x       nan     nan
006     0          3.00   19.00
007     0          3.00   13.00
008     0          9.00   13.00
009     0          9.00   19.00
010     0          3.00   19.00
011     0   o      8.00   18.00
012     0          8.00   14.00
013     0          4.00   14.00
014     0          4.00   18.00
015     0          8.00   18.00
016     0   x       nan     nan
017     0          6.00   17.00
018     0          5.00   15.00
019     0          7.00   15.00
020     0  ___     6.00   17.00
021     1   o     12.00   18.00
022     1         12.00   12.00
023     1         20.00   12.00
024     1         20.00   10.00
025     1         10.00   10.00
026     1         10.00   20.00
027     1         20.00   20.00
028     1         20.00   18.00
029     1         12.00   18.00
030     1   o     25.00   24.00
031     1         25.00   14.00
032     1         15.00   14.00
033     1         15.00   16.00
034     1         23.00   16.00
035     1         23.00   22.00
036     1         15.00   22.00
037     1         15.00   24.00
038     1  ___    25.00   24.00
039     2   o     14.00   20.00
040     2         10.00   20.00
041     2         15.00   28.00
042     2         14.00   20.00

 

s2.IFT

array([[ 0,  0, 11],
       [ 0, 11, 21],
       [ 1, 21, 30],
       [ 1, 30, 39],
       [ 2, 39, 43]])

 

Arcpy geometry representation

 

This is the dissection of the first polygon down to its elemental parts and the arcpy class methods and properties that can be accessed through the standard interface.

 

arcpy.da.SearchCursor

cur = arcpy.da.SearchCursor(in_fc2, 'SHAPE@', spatial_reference=SR)
dir(cur)
[[... snip ..., '_as_narray', '_dtype', 'fields', 'next', 'reset']

arcpy.Polygon

p0  # ---- the first polygon
<Polygon object at 0x1e5284a3320[0x1e5284c1e18]>
dir(p0)
['JSON', 'WKB', 'WKT', … '_fromGeoJson', …  'angleAndDistanceTo', 'area', 'boundary', 'buffer', 'centroid',
 'clip', 'contains', 'convexHull', 'crosses', 'cut', 'densify', 'difference', 'disjoint', 'distanceTo', 'equals', 'extent',
 'firstPoint', 'generalize', 'getArea', 'getGeohash', 'getLength', 'getPart', 'hullRectangle', 'intersect', 'isMultipart',
 'labelPoint', 'lastPoint', 'length', 'length3D', 'measureOnLine', 'overlaps', 'partCount', 'pointCount',
 'pointFromAngleAndDistance', 'positionAlongLine', 'projectAs', 'queryPointAndDistance', 'segmentAlongLine',
 'snapToLine', 'spatialReference', 'symmetricDifference', 'touches', 'trueCentroid', 'type', 'union', 'within']

arcpy.Array

p0[0]  # ---- first polygon's first part... aka, an array of point objects
<Array [<Point (300010.0, 5000020.0, #, #)>, … snip, <Point (300003.0, 5000019.0, #, #)>]>
dir(p0[0])
[ … snip ..., 'add', 'append', 'clone', 'count', 'extend', 'getObject', 'insert', 'next', 'remove', 'removeAll', 'replace', 'reset']
arcpy.Point
dir(p0[0][0])  # ---- the point
['ID', 'M', 'X', 'Y', 'Z', … snip …, 'clone', 'contains', 'crosses', 'disjoint', 'equals', 'overlaps', 'touches', 'within']

...the coordinates

p0[0][0].X  # ----  finally, the X coordinate of the first point of the first array in the first shape of the
                            first polygon... got it?
300010.0

Geo class vs ndarray

This is a summary of the methods and properties that I have added to the Geo array
gs = set(dir(g))
ss = set(dir(g.base))
sorted(gs.difference(ss))
['AOI_extent', 'AOI_rectangle', 'FT', 'IDs', 'IFT', 'Info', 'K', 'N', 'X', 'XY', 'Y', 'Z', '__dict__', '__module__',
 'angles', 'areas', 'bits', 'centers', 'centroids', 'convex_hulls', 'densify_by_distance', 'extent_rectangles', 'extents',
 'fill_holes', 'get', 'holes_to_shape', 'is_convex', 'is_multipart', 'lengths', 'maxs', 'means', 'min_area_rect', 'mins',
 'move', 'od_pairs', 'outer_rings', 'part_cnt', 'parts', 'pnt_info', 'pnts', 'polylines_to_polygons', 'polys_to_segments',
 'pull', 'rotate', 'shapes', 'split', 'translate', 'unique_pnts']

Continued at...

Options for different representations of arcpy geometry arrays are there.

I will continue the development of the Geo class based on numpy's ndarray in my GitHub at...

 

npGeo ... a geometry class

Dan_Patterson

The Py... Links II

Posted by Dan_Patterson Champion May 23, 2019

The Python links II.

 

The original Py Links was getting a bit packed, and with the demise of python 2 on the horizon, I thought I would make a clean break and focus on python 3.x and its kindred as they ship with ArcGIS Pro.  This link will be a bit thin for a while, but it will fill I am sure.

 

This is a new listing of things pythonic as they relate to GIS analysis.  It is organized largely by theme and/or package.

 

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

New

 

Which packages will require python 3? and which ones will support 2.7 and 3.x after December ... 2019-10-08

 

Python 3.6.9 …. security fixes only... last of the 3.6 line prep for ArcGIS Pro 2.5... 2019-09-25

 

Machine Learning Algorithms …. for the beyond mapping types.... 2019-08-27

 

Python readiness check... 2019-08-27

 

The Walrus operator in Python 3.8 ... (PEP 572)  : 2019-07-30 

Cool operator, cause of the rift with the BDFL and the rest, perhaps an homage to the Beatles

 

Speedup your algorithms... Pytorch, Numba, Parallelization, Dask   :  2019-07-20

An excellent 4 part series on improving algorithm speed.  The approaches are different and often complementary.

 

Spyder 3.3.5 and ArcGIS 1.6.2 available   : 2019-07-16 

Update spyder and ArcGIS through conda if installed that way.

As part of the install https://json5.org/Json5 (json for humans) is also installed

 

Python 3.7.4 final release available           : 2019-07-08

Live on the edge... only for Pro though

 

 

 

: ---------------------------------------------------------------------------------------------------------------------------- :

Python

News

JSON5 ... json for humans  2019-07-16

Kite - AI Autocomplete for python coding 2019-06-21  May be coming to Spyder

keras feature extraction on large datasets with deep learning  2019-05-28

New features planned for python 4.0  2019-05-23

Yes, python 4.0.  Still on 2.x???

 

Version and changelog

Changelog — Python 3.7.2 documentation :

The most current version will be listed, you can descend the tree to find out when particular aspects were implemented.

What’s New in Python — Python 3.8.0a4 documentation 

Same as above, the most current or development version is listed first with previous versions descending.

IDEs for scripting

The ones listed here are available for installation in either your base ArcGIS Pro installation path or a clone of it, depending on whether you are the master of your universe at home/work.

Clone... ArcGIS Pro ... for non administrators

Spyder

Spyder.... for coding with Python :  This link provides a visual guide to some of the functionality within Spyder

Spyder docs...  official documentation.

Spyder on GitHub...  A good place to follow changes.

Qt Console

Qt Console... ships with anaconda

Jupyter Lab

Github page, another Jupyter project

Jupyter Notebook

The notebook homepage

Optimize Jupyter Notebooks...               Cool article, some tips apply to spyder

: ---------------------------------------------------------------------------------------------------------------------------- :

Numpy

NumPy docs...

Release notes...

NumPy on GitHub...

: ---------------------------------------------------------------------------------------------------------------------------- :

SciPy

scipy docs...

SciPy on GitHub ...

: ---------------------------------------------------------------------------------------------------------------------------- :

Pandas

pandas docs...

Pandas on GitHub ...

: ---------------------------------------------------------------------------------------------------------------------------- :

ArcGIS Pro

System requirements ...

ArcGIS Pro 2.4 release notes...

: ---------------------------------------------------------------------------------------------------------------------------- :

Old News

is good news

 

Kite - AI Autocomplete for python coding    : 2019-06-21  May be coming to Spyder

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

Export a single shape from a featureclass.  

Looks the same ... right?

It isn't, necessarily.

 

Here is an example:

 

A single multipart shape.

The point order would be the same, wouldn't it?

 

Not sure any more... some different rule set seems to be applied when a shape leaves its 'nest' and ventures out into the world on its own.

 

I can understand point order being influenced by the construction method and its association with other shapes.  In the case of the shape 'in-group', the outer part is ordered from the top right(10, 20) and clockwise.  The second, smaller part is the listed second (line 13).

 

When Shape 1, ventures out on its own, the inner part is now first (line 2), almost like freedom provide a new rule-set during its transition

 

Shape 1... its in-group pointsShape 1... out on its own
array([
array([[10., 20.],
       [10., 10.],
       [ 0., 10.],
       [ 0., 20.],
       [10., 20.],
       [nan, nan],
       [ 2., 18.],
       [ 2., 12.],
       [ 8., 12.],
       [ 8., 18.],
       [ 2., 18.]]),
array([[ 7., 17.],
       [ 7., 13.],
       [ 3., 13.],
       [ 3., 17.],
       [ 7., 17.],
       [nan, nan],
       [ 5., 16.],
       [ 4., 14.],
       [ 6., 14.],
       [ 5., 16.]])
], dtype=object)
array([
array([[ 7., 17.],
       [ 7., 13.],
       [ 3., 13.],
       [ 3., 17.],
       [ 7., 17.],
       [nan, nan],
       [ 5., 16.],
       [ 4., 14.],
       [ 6., 14.],
       [ 5., 16.]]),
array([[10., 20.],
       [10., 10.],
       [ 0., 10.],
       [ 0., 20.],
       [10., 20.],
       [nan, nan],
       [ 2., 18.],
       [ 2., 12.],
       [ 8., 12.],
       [ 8., 18.],
       [ 2., 18.]])
], dtype=object)

So, there is probably a logical explanation for what is seen, BUT, if you were relying on the initial point ordering when the shape lived at home, to apply when it moved out on its own, then you would be in for a surprise.

 

Conclusions:

 

  • Don't rely on what you see.
  • Examine what you have.
  • Be prepared to 'deal'
  • Not all shapes behave the same way when they move out.
  • Any rules that you come up with probably have a corner-case
  • At least it didn't take anything that didn't belong to it, when it moved out