2019

May 2019

# Geometry : Part 8

Posted by Dan_Patterson 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.IFTarray([[ 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

Posted by Dan_Patterson May 23, 2019

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

Python 3.9 beta... lots of new stuff ... 2020-04-08  perhaps a long wait in Arc*** world

2020-04-22  A great introduction for handling paths

Python cheat sheet ... 2020-04-02  An really good cheat sheet!  One of the best I have seen

Arcgis 1.8 release notes ...  2020-03-09

Panda 1.0.0 released ... 2020-30-01

NumPy, the absolute beginners guide 2020-01-21

Panda 1.0.0... what's new ... 2020-01-15

PYTHON 2.7 RETIRES ... 2019-12-31

Why is the migration to python 3 taking so long  An interesting history as well  2019-11-15

Guido is heading into retirement Pythonistas.... Guido is heading into retirement 2019-11-05

NumPy 1.17.4 released  Conda install in your ArcGIS Pro distribution or clone works fine.  2019-11-03

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

https://blog.datasciencedojo.com/machine-learning-algorithms/

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

# Python

## News

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

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

New features planned for python 4.0  2019-05-23

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

## Version and changelog

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

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

Cool article, some tips apply to spyder

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

# Numpy

NumPy docs...

Release notes...

NumPy on GitHub...

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

# SciPy

scipy docs...

SciPy on GitHub ...

Scipy .... Nature article a great read

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

# 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

for the beyond mapping types.... 2019-08-27

2019-08-27

: 2019-07-30

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

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

: 2019-07-08

Live on the edge... only for Pro though

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

Posted by Dan_Patterson May 13, 2019

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)shapesarray([(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

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

Posted by Dan_Patterson May 8, 2019

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

#### Filter Blog

By date: By tag: