Skip navigation
All People > Dan_Patterson > Py... blog
1 2 3 Previous Next

Py... blog

120 posts

Cool

 

It is nice that you can view geometry in ArcGIS Pro.

 

Ditto for notebooks in your browser.

 

But I really hate cranking up a new featureclass when I am working on a geometry exploration, when all I want to see is what the numbers actually look like.

 

I stumbled on this when I was working on my npgeom package, which uses an alternate geometry constructor than is used in the arc* line of products.  It also deconstructs geometry using arcpy data access cursors and/or  FeatureClassToNumPyArray.   (Thas is another story though)

 

In short, I was doing my thing and got 3 polygons objects from a featureclass.

Lines 78-80... blah blah, stuff comes out

Line 81 ... wanted to make sure they were polygon objects.... indeedy they were

Line 82, 83, 84 (right side of graphic)

   cool!  perfect polygons

       [82] A multipart polygon with 2 holes in the first part and 1 in the second

       [83] Another multipart, 3 holes in one, none in the other

       [84] The triangle, paying homage to a simpler time when geometry deconstruction was easier

 

 

I won't go into details, but my python setup is noted at the bottom of the callout in the graphic.

 

More details when I explore more.  Going to replicate this as a *.ipynb for use in the browser AND with Spyder as well.

Hopefully other python IDEs support these as well.  Some use qt and mpl in their graphics display arsenal.

 

Addendum

I did forget to mention that you can save the contents of the qt console in spyder to a couple of formats for posterity.

I save a little sample in the attached zip file.. 

  • unzip it to a location (eg. c:\temp
  • double-click on the npg_01.html file and it should load in your browser

Of course you can edit the html file to fix any stuff that you want.

 

You can save to svg format as well (see attached)

 

You people had better get back to work

 Free basic functionality.  

 

Once again, functionality that is normally restricted to the Standard or Advanced ArcGIS Pro license. 

 

 

Previously

 - Feature to Point 

 - Feature Extent to Poly Features

 - Free Tools : Frequency and Statistics

 - Free Tools : Convex Hulls

 

Help topics

 - Feature Envelope to Polygon

 - Frequency

 -  Convex hull in Minimum bounding geometry

 -  Feature to Point

 -  Polygons to line

 -  Circle in Minimum bounding geometry

 

Implementation

The implementation here is largely based on Welzl's algorithm. 

Another container.  I did the rectangular ones, now the circle, then the ellipse.

Still, numpy, python and arcpy all play nice.  Check out the code.  All seven tools are implemented in one toolbox and controlled by one script.

 

Output examples

A circle is implemented as an ngon with 180 sides (2 degree increments) as polygons.  

I could have put in a bunch of options in the tool to select the density and output type, but students (and others) wouldn't learn to 'tweak' code.  You can always do your own toolbox or use esri's

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

 

WARNING

I have code that checks the validity of the file paths.  If you input or output paths contain spaces or other flotsam, then the tool will not produce any results.  Why?  Too many questions where file paths are the problem.  I won't 'enable' the current practice

 

Download

You can copy the contents of the folder on my GitHub pages.  No fancy install, just create a folder, throw the stuff in, load the toolbox and give it a try.

 

Got any geometry related or analysis tools you need implemented? let me know

 

Free_Tools

 Free basic functionality.  

 

Once again, functionality that is normally restricted to the Standard or Advanced ArcGIS Pro license. 

 

 

Previously

 - Feature to Point 

 - Feature Extent to Poly Features

 - Free Tools : Frequency and Statistics

 - Free Tools : Convex Hulls

 

Help topics

 - Feature Envelope to Polygon

 - Frequency

 -  Convex hull in Minimum bounding geometry

 - Feature to Point

 - Polygons to line

 

Implementation

The implementation here is a combination of a couple of things. 

  • conversion of Polygons to Polylines
  • poly* features to segments

It doesn't do the overlap thing... you will have to pay the big $$$ for that (for now).

 

Output examples

I tossed in a bounding container so you could see the points. 

  • Specify the input featureclass (polygon, polyline or multipoint),
  • the output featureclass,
  • select feature to point from the tool selection and

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

 

 

 

 

First... To the left.

Not too exciting, but it is there.

 

Second... To the right.

Pretty.

 

 

 

 

WARNING

I have code that checks the validity of the file paths.  If you input or output paths contain spaces or other flotsam, then the tool will not produce any results.  Why?  Too many questions where file paths are the problem.  I won't 'enable' the current practice

 

Download

You can copy the contents of the folder on my GitHub pages.  No fancy install, just create a folder, throw the stuff in, load the toolbox and give it a try.

 

Got any geometry related or analysis tools you need implemented? let me know

 

 

 

Free_Tools

 Free basic functionality.  

 

Once again, functionality that is normally restricted to the Standard or Advanced ArcGIS Pro license. 

 

 

Previously

 - Feature Extent to Poly Features

 - Free Tools : Frequency and Statistics

 - Free Tools : Convex Hulls

 

Help topics

 - Feature Envelope to Polygon

 - Frequency

 -  Convex hull in Minimum bounding geometry

 - Feature to Point

 

Implementation

The implementation here is basically the default. Using the average or weighted average of the points making up the feature.  Soooo basically some kind of average, like you can get the points yourself using ...FeatureClassToNumPyArray… with the 'explode_to_points' option.  Voila! the points and their coordinates, ergo, the averages in some form (weighted or simple)

 

Output example

I tossed in a bounding container so you could see the points. 

  • Specify the input featureclass (polygon, polyline or multipoint),
  • the output featureclass,
  • select feature to point from the tool selection and

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

 

 

 

The variant shown here is showing the full shape as input.  A tweak, enables you to separate out multipart shapes if you want, or more simply use the ...Multipart to Singlepart... tool if you have a crushing need to carry over the attributes... It will save you a Join and arcpy has to do something occasionally.

 

The conversion uses the numpy based Geo class that I describe in the 8 part series on geometry. 

You could add the geometry attributes to the result if needed ...Add Geometry Attributes … 

 

WARNING

I have code that checks the validity of the file paths.  If you input or output paths contain spaces or other flotsam, then the tool will not produce any results.  Why?  Too many questions where file paths are the problem.  I won't 'enable' the current practice

 

Download

You can copy the contents of the folder on my GitHub pages.  No fancy install, just create a folder, throw the stuff in, load the toolbox and give it a try.

 

Free_Tools

 

Got any geometry related or analysis tools you need implemented? let me know

 Free basic functionality.  

 

Another Free missive that uses numpy and arcpy to produce functionality that is normally restricted to the Standard or Advanced ArcGIS Pro license. 

 

 

Previously

 - Feature Extent to Poly Features

 - Free Tools : Frequency and Statistics

 

Help topics

 - Feature Envelope to Polygon

 - Frequency

 -  -Convex hull in Minimum bounding geometry

 

Spatial containers

The extent poly* features done previously, is one of the standard containers.  The convex hull is the most widely used and the easiest to implement, not because of the simplicity but because of the availability of standard algorithms.  Many packages use modules from the ... qhull … package.

 

Normally containers only make sense if you are using projected coordinates or can perform geodesic densification. 

 

Output example

Pretty well sums it up. 

  • Specify the input featureclass (polygon, polyline or multipoint),
  • the output featureclass,
  • select convex hulls from the tool selection

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

 

 

 

The conversion uses the numpy based Geo class that I describe in the 8 part series on geometry. 

 

I could add the original attributes to the result (either within the toolset or after) or I could also 

Add Geometry Attributes ...

 

The full call to the tool, or the equivalent bits that I need.

A spatial or attribute join would be another alternative if you need attributes as well.

 

If you have a preference let me know.

The results are derived quickly and there is an optimization if the number of points making up the shape exceed about 50.  This was a qualitative estimate of the cross-over point between implementing a python solution versus a C compiled solution from qhull.

 

 

WARNING

I have code that checks the validity of the file paths.  If you input or output paths contain spaces or other flotsam, then the tool will not produce any results.  Why?  Too many questions where file paths are the problem.  I won't 'enable' the current practice

 

Download

You can copy the contents of the folder on my GitHub pages.  No fancy install, just create a folder, through the stuff in, load the toolbox and give it a try.

 

Free_Tools

A dministrator privileges … or you know the IT peeps … or you have created a cloned environment.

Pick one.

 

My installation path :     C:\arc_pro   ….. everything beyond this point is the same

Your installation path :   C:\...........    ….. got it?

 

Table of contents

 

Download and install tips

 

1  Follow the help topics:

ArcGIS Pro 2.4 system requirements—ArcGIS Pro | ArcGIS Desktop 

Download, install, and authorize—ArcGIS Pro | ArcGIS Desktop 

 

2  Go back to step 1.

Really, it is good and should be read, especially the part about your computer being able to run the software

 

3  My Esri, My Organization, Downloads

If it is there, it will look like the following:

 

 

4  Installation steps for retentives

Now, don't hit the Run option!  It is tempting, but there is Save and Save As.  Save As will be used.

To prepare for this, you should have done the following (not!, I am guessing)

  •  Make a folder to download your software ....
    • C:\users\you\whatever\downloads ... is no good, just because
    • C:\Computer\ArcGISPro_24 is good... simple, obvious and you own it
    • Download the *.exe to that folder using Save As
  • Right-click on the *.exe file and run it as administrator, specifying the above folder as the destination
  • Do the same for the *.msi file and you will automagically get a bunch of stuff in that folder AFTER the installation is complete... just follow that, but your folder should look like the following

Where step 1 is the main installation folder you created and downloaded the *.exe (2), when you run the *.exe, you will get the folder in step 3, and run the *.msi and you get the rest of the stuff.

 

Why do I do this? 

Because if things go really really bad, you will know where the ArcGISPro.msi file is, so when you have to do a complete uninstall, you can reinstall within a minute. 

Simple... no remembering or letting Microsoft Parent decide where things should go

 

What I did next

I do the conda thing... some legacy but relevant reading

ArcGIS Pro 2... Creating Desktop Shortcuts 

Spyder.... for coding with Python 

ArcGIS PRO  .... your conda environments and script editor 

 

Crank up conda through whatever means to run ...proenv.bat which sets everything up.  What is show below is what happens when I created a shortcut (Dolly) and messed around with the python ide so it isn't as dark and gloomy as yours will be.

 

 

I needed the following to do the programming I need and I did it in the following order.

 

1  Update numpy

(arcgispro-py3) C:\arc_pro\bin\Python\envs\arcgispro-py3> conda update numpy

 

2 Downgrade sphinx to 1.8.5  (needed IF you document your scripts, otherwise the documentation will look horrible)

(arcgispro-py3) C:\arc_pro\bin\Python\envs\arcgispro-py3> conda install sphinx==1.8.5

 

3  installed sphinx_rtd_theme   Getting Started with Sphinx — Read the Docs 3.5.3 documentation 

    You can skip this step if you don't do documentation or produce reports, or use Markdown or reStructured Text (or know what I am talking about )

(arcgispro-py3) C:\arc_pro\bin\Python\envs\arcgispro-py3>conda install sphinx_rtd_theme --no-pin

 

4  Install spyder

(arcgispro-py3) C:\arc_pro\bin\Python\envs\arcgispro-py3>conda install spyder

 

Tips

Never, never install without doing a test run first!

 

      (arcgispro-py3) ….snip …. >conda install some_package --dry-run 

 

Then examine what it is going to do.  Sometimes, nothing 'bad' will happen, but you should at least make a copy what you are about to install.  If things go bad, you can roll back through the 'revisions' to a previous state.

 

Revision History from this install

 

(arcgispro-py3) C:\arc_pro\bin\Python\envs\arcgispro-py3>conda list --revisions
2019-06-27 20:36:30  (rev 0)  Fresh install of ArcGIS Pro 2.4
    +arcgis-1.6.1 (esri)
    +arcgispro-2.4 (esri)
    ... huge snip ....
    +zeromq-4.3.1
    +zlib-1.2.11

2019-06-27 20:40:06  (rev 1)    The numpy upgrade
     ca-certificates  {2019.1.23 -> 2019.5.15}
     certifi  {2019.3.9 -> 2019.6.16}
     cffi  {1.12.2 -> 1.12.3}
     .... snip ...              
     numpy  {1.16.2 -> 1.16.4}
     numpy-base  {1.16.2 -> 1.16.4}
     .... snip .... 
    +pywin32-223
    +zipp-0.5.1
2019-06-27 20:47:46  (rev 2)   And So On.
    +alabaster-0.7.12
 .... snip .... 
2019-06-27 22:07:03  (rev 4)  And finally
    +sphinx_rtd_theme-0.4.3

Now if anything goes wrong, (Assuming I want to go back to revision 1)

(arcgispro-py3) C:\arc_pro\bin\Python\envs\arcgispro-py3>conda install --revision 1  (change 1 to your revision)


A little conda in spyder anyone?

Note:
There are load of IPython line and cell magics that can be used with Spyder.
Summary of magic functions (from %lsmagic):

Available line magics:
  %aimport  %alias  %alias_magic  %autoawait  %autocall  %automagic  %autoreload
  %autosave  %bookmark  %cd  %clear  %cls  %colors  %conda  %config  %connect_info
  %copy  %ddir  %debug  %dhist  %dirs  %doctest_mode  %echo  %ed  %edit  %env  %gui
  %hist  %history  %killbgscripts  %ldir  %less  %load  %load_ext  %loadpy  %logoff
  %logon  %logstart  %logstate  %logstop  %ls  %lsmagic  %macro  %magic  %matplotlib
  %mkdir  %more  %notebook  %page  %pastebin  %pdb  %pdef  %pdoc  %pfile  %pinfo
  %pinfo2  %pip  %popd  %pprint  %precision  %prun  %psearch  %psource  %pushd
  %pwd  %pycat  %pylab  %qtconsole  %quickref  %recall  %rehashx  %reload_ext  %ren
  %rep  %rerun  %reset  %reset_selective  %rmdir  %run  %save  %sc  %set_env  %store
  %sx  %system  %tb  %time  %timeit  %unalias  %unload_ext  %varexp  %who  %who_ls
  %whos  %xdel  %xmode

Available cell magics:

  %%%%HTML  %%SVG  %%bash  %%capture  %%cmd  %%debug  %%file  %%html  %%javascript
  %%js  %%latex  %%markdown  %%perl  %%prun  %%pypy  %%python  %%python2  %%python3
  %%ruby  %%script  %%sh  %%svg  %%sx  %%system  %%time  %%timeit  %%writefile

Certainly enough command line stuff to reminisce about the days of 40 character displays and green crts

 

Good luck

F ree basic functionality.  

 

I use Free as both a verb and an adjective

The examples I will be providing should be Free(d) and should not require an Advanced or Standard License... they should be Free with a basic license.

 

 

Previously

 - Free Tools : Frequency and Statistics

 

Help topic

Feature Envelope to Polygon

 

Spatial containers

Spatial containers can be used as substitutes for the original spatial pattern.  These would include ...

  • the shape itself,
  • convex hull,
  • concave hull
  • minimum area bounding (or perimeter)
    • circle
    • ellipse
    • rectangle
  • extent poly* features

 

The last of these are axis-aligned shapes bounded by the coordinates of the minimum of the lower left corner to upper right corner, not of the existing coordinates, but by the L(eft), B(ottom), R(ight) and T(op)… LBRT. 

Normally containers only make sense if you are using projected coordinates or can perform geodesic densification. 

 

Output example

Pretty well sums it up. 

  • Specify the input featureclass (polygon, polyline or multipoint),
  • the output type (polygon or polyline) and
  • the output featureclass.

The conversion uses the numpy based Geo class that I describe in the 8 part series on geometry. 

 

The conversion also handles shapes that contain curves, for simple geometry.  To do this, you will find that circles and ellipses or sliced sphericals contain 2 points... the identical start and end and no other points.  At least in a file geodatabase.  The quick solution is to densify the arc based on the ANGLE option in arcpy.densify.  I choose between 1, 2 or 5 degree densification depending on how fine I want the resultant n-gon to reflect the original  curve.

 

I could add the original attributes to the result (either within the toolset or after) or I could also 

Add Geometry Attributes ...

The full call to the tool, or the equivalent bits that I need.

A spatial or attribute join would be another alternative if you need attributes as well.

 

If you have a preference let me know.

The script is embedded for now.  I will release the script and final version when I have addressed all issues.

 

 

WARNING

I have code that checks the validity of the file paths.  If you input or output paths contain spaces or other flotsam, then the tool will not produce any results.  Why?  Too many questions where file paths are the problem.  I won't 'enable' the current practice

F ree basic functionality.  

 

I use Free as both a verb and an adjective

The examples I will be providing should be Free(d) and should not require an Advanced or Standard License... they should be Free with a basic license.

History

Frequency .... the help topic

It started a long time ago with Split Layer By Attributes.

Advanced license for so long, when really... the means to provide the analysis were already there.  Eventually, it was Free(d).

 

So this will be a series on how to perform functionality using the tools that are already provided for you within ArcGIS Pro and supplemental Python Packages which are already installed.

 

Maybe some educators out there will take application development directed towards analysis  as a serious area of GIS.

 

Tools to demonstrate

This blog series will NOT be about making maps or stuff related to expedite map making.  So if you want to develop tools for map making, then look elsewhere. 

I will begin with tools in the Analysis and Data Management Tools, for example:

 

  • Analysis Tools 
    •  Thiessen PolygonsThiessen Polygons in ArcGIS Pro  (Actually, I have done several incarnations of spatial triangulation and allocation already... nobody notices...)
  •  Data Management Tools
    •  Frequency … Frequency tool  (Ditto, One exists on the code sharing site but this is a beefed up version)

 

Example of Frequency and Statistics

If you had 4 counties with 4 towns of various size classes and you were interested in studying age dynamics versus population size, you might want to start with that basic premise.

  • Classify/group/categorize the towns by population size
  • Produce a unique classification scheme or use 2 or more existing fields to group your data for counting (aka, frequency).
  • Now, that the data are grouped, you can now determine some basic statistics for those classes.

 

Output example

See the table as an example...  A table is as good as a map and you don't get distracted by all those colors.

 

 

So, row 1, County A, Town_class A_.... 1195 people, some stats... Age_min will not be below 18... pretend, survey privacy.

 

A little Chi-Square, perhaps a Moran's (if we wanted to map) and off you go.  

Sadly this is in an Advanced license.  The whole production of the classification scheme, the frequency determination and the sorting is basically 1 line of code.  The stats stuff? python and numpy.  Null values? handled with ease... besides you should never have nulls anyway.

 

Requirements

So download and try it out on your data.

I would be interested in use cases to see what might be added.

The only restrictions...

  • ArcGIS Pro 2.4+ (might work on lower versions, but I no longer have them, and these are free anyway)
  • Locally stored data in a file geodatabase, like gdb tables or featureclass tables
  • If you have excel or csv files... do the work and make them gdb tables (there are tools for that, check ArcToolbox or my blog)
  • I embedded the script into the toolbox, let me know if it doesn't work embedded.  When I finish with suggestions for additions, I will Free the script so you can make your own modifications or additions

 

Have fun.

 

Up Next

Extent to Polygons  Converting feature extent(s) to polygons is a relatively easy task 

Featureclass Properties

 

  • Feature class info
  • Field info
  • Geometry info
  • Geometry decomposition

 

The code... I will update here.

 

 

Information functions for featureclasses

A quick solution to obtain geometry information.  Other info can be added.

fc_info(in_fc, prn=True)

FeatureClass:
   C:/Arc_projects/CoordGeom/CoordGeom.gdb/Shape2
shapeFieldName  OIDFieldName  shapeType spatialReference
Shape           OBJECTID      Polygon   NAD_1983_CSRS_MTM_9

Those pesky fields.  Don't want to open up ArcGIS Pro to find out?  Already got Spyder (or your Python IDE) open? Use this.

fld_info(in_fc, prn=True)

FeatureClass:
   C:/Arc_projects/CoordGeom/CoordGeom.gdb/Shape2
Name          Type         Length Nullable  Required 
OBJECTID      OID               4 False     True     
Shape         Geometry          0 True      True     
Shape_Length  Double            8 True      True     
Shape_Area    Double            8 True      True     
CENTROID_X    Double            8 True      False    
CENTROID_Y    Double            8 True      False    
INSIDE_X      Double            8 True      False    
INSIDE_Y      Double            8 True      False    

How many shapes? Are they all singlepart? How many points?  Which points connect to make what shape?

fc_geom_info(in_fc, SR=None, prn=True, start=0, num=10)

Featureclass:
    C:/Arc_projects/CoordGeom/CoordGeom.gdb/Shape2
   Shape    Parts   Points From_pnt   To_pnt
       1        2       21        0       21
       2        2       18       21       39
       3        1        4       39       43

Same as above... but where do the parts stop and start? Everyone loves stats! Right?

fc_composition(in_fc, SR=None, prn=True, start=0, end=50)

C:/Arc_projects/CoordGeom/CoordGeom.gdb/Shape2
Shapes :   3
Parts  :   5
  max  :   2
Points :   43
  min  :   4
  median : 9
  max  :   11
     IDs     Part   Points From_pnt   To_pnt
       1        0       11        0       11
       1        1       10       11       21
       2        0        9       21       30
       2        1        9       30       39
       3        0        4       39       43

 

So simple, so...  _common give it a try

More offerings in … npGeo ... numpy geometry 

 

Oh yes... Nice documentation too...

 

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 

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.2 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

 

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  
Geometry

----

As before, the inputs

----The polygons that I will be using are shown to the right.
  1. A square, 5 points, first and last duplicates
  2. A lake on an island in a lake...
  3. A multipart with a two different shaped donut holes
  4. The letter 'C'
  5. The letter 'V'
Each part is labelled at the labelpoint rather than the centroid, hence each part gets labelled.

 

 

 

 

 

---- An Alternate Geometry Reconstructor ----

---- Poly* to arrays revisited ----

 

These two scripts help blow poly features into their constituent parts.  Initially a null point (np.nan, np.nan) is used to separate the poly features.  The location of these insertions is retained and a fr_to index list is maintained so that the resultant 2D points array can be reconstructed or used in other calculations.

 

The outputs are: a_2d, ids, fr_to, id_fr_to which represent the 2D array, the id values of where one poly feature ends, a from-to list is constructed and a final list with the ids prepended is also included.

def _pp_(poly):
    """See poly_pnts for details.  This is for single poly feature conversion
    requires
        null_pnt = (np.nan, np.nan)
    """

    sub = []
#    for i, arr in enumerate(poly):
    for arr in poly:
        pnts = [[pt.X, pt.Y] if pt else null_pnt for pt in arr]
        sub.append(np.asarray(pnts)) #append(pnts)
    return np.asarray(sub)


def poly_pnts(in_fc, as_obj=False):
    """Poly features to points while maintaining the location of the points in
    the input features. null points are replaced with their numpy equivalent.

    Parameters
    ----------
    in_fc : text
        Full path to the featureclass.
    as_obj : boolean
        True returns a list of arrays of the shapes.  The array types may be
        mixed.  False, returns a 2D array of points and an array of indices.

    Returns
    -------
    False, returns a 2D array of points and the location to split
    those points should you need to reconstruct the poly feature.  A `nan`
    point separates the parts of the features within the poly feature, so split
    first. If `as_obj` is True, a list of arrays is returned. 

    Notes
    -----
    Use `polys = fc_shapes(in_fc)` to obtain the poly features separately.
    """

    SR = getSR(in_fc)   
    polys = []
    with arcpy.da.SearchCursor(in_fc, 'SHAPE@', spatial_reference=SR) as cur:
        polys = [row[0] for row in cur]
    id_too = []
    a_2d = []
    for i, p in enumerate(polys):
        r = _pp_(p)                              # calls to _pp_
        id_too.extend([(i, len(k)) for k in r])
        a_2d.extend([j for i in r for j in i])
    a_2d = np.asarray(a_2d)
    id_too = np.array(id_too)
    ids = id_too[:, 0]
    too = np.cumsum(id_too[:, 1])
    frum = np.concatenate(([0], too))
    fr_to = np.array(list(zip(frum, too)))
    id_fr_to = np.array(list(zip(ids, frum, too)))
    if as_obj:
        a_obj = [a_2d[f:t] for f,t in fr_to]  # alternate constructor
        return a_obj, ids, fr_to, id_fr_to
    return a_2d, ids, fr_to, id_fr_to
           
   

 

 

A sample run

  • Line 1 reads the featureclass, produces the 2D point array and the various indices
  • In Line 2 and 3, I usually calculate the mean of the point cloud accounting for the null points separating the individual poly features.
  • A separate object array can be constructed if you have a need to work with the individual poly features at once, otherwise, you can recreate this using the a array and the fr_to indices as shown in line 4. 
  • Line 5 reconstructs the original polygons so that they can be moved back into ArcGIS Pro  ... more later

Normally you would 'do stuff' between line 4 and 5 and send back the resultant array, but this is just for illustration

 

a, ids, fr_to, id_fr_to = poly_pnts(in_fc)
m = np.nanmin(a, axis=0)
a_s =- m
a1 = np.asarray([a[f:t] for f,t in fr_to])
p_arr = [_arr_poly_(i, SR) for i in a1]  # ** to reverse np.concatenate(a1)
frmt = """
Polygon ids:   {}
From-to pairs:
{}
Id_from_to array
{}
"""

print(dedent(frmt).format(ids, fr_to, id_fr_to))

Polygon ids:   [0 1 1 2 2 3 4]
From-to pairs:
[[ 0  5]
[ 5 16]
[16 26]
[26 36]
[36 48]
[48 57]
[57 65]]
Id_from_to array
[[ 0  0  5]
[ 1  5 16]
[ 1 16 26]
[ 2 26 36]
[ 2 36 48]
[ 3 48 57]
[ 4 57 65]]

Line 15 shows that poly features 1 and 2 have 2 parts.

Lines 16-23 are the from-to pairs of points needed to reconstruct the arrays to make polygons.

Line 24- just combines the two previous lists.

Reconstructing the arrays to poly features

The code below does this.  A helper function, then the code block that uses a search cursor to reassemble the array to something that arcpy can use.  Sadly you have to go from a numpy array, then create points, which are placed in an arcpy Array and from there the arcpy.Arrays are assembled to form polygon or polyline features.  And Finally!!! A single-part featureclass is created, then a multipart featureclass, like the original data that went into the whole process. 

 

Complicated? Not really, the real bottleneck is the cursors.  You need them going in (regardless of the shroud placed around them) and going out. (ie __geo_interface__,  _as_narray,  FeatureClassToNumPyArray.

I should point out here that Numpyarraytofeatureclass works with simple geometry to get poly features back, but who works with simple features all the time.

 

def _arr_poly_(pnts, SR):
    """Single array to polygon, array can be multipart with or without interior
    rings.  Outer rings are ordered clockwise, inner rings (holes) are ordered
    counterclockwise.  For polylines, there is no concept of order
    Splitting is modelled after _nan_split_(arr)
    """

    subs = []
    s = np.isnan(pnts[:, 0])
    if np.any(s):
        w = np.where(s)[0]
        ss = np.split(pnts, w)
        subs = [ss[0]]
        subs.extend(i[1:] for i in ss[1:])
    else:
        subs.append(pnts)
    aa = []
    for sub in subs:
        aa.append([arcpy.Point(*pairs) for pairs in sub])
    poly = arcpy.Polygon(arcpy.Array(aa), SR)
    return poly

      
def arr_poly_fc(a, p_type='POLYGON', gdb=None, fname=None, sr=None, ids=None):
    """Reform poly features from the list of arrays created by poly_pnts

    Parameters
    ----------
    a : array or list of arrays
        Some can be object arrays, normally created by ``pnts_arr``
    p_type : string
        Uppercase geometry type
    gdb : text
        Geodatabase name
    fname : text
        Featureclass name
    sr : spatial reference
        name or object
    ids : list/array
        Identifies which feature each input belongs to.  This enables one to
        account for multipart shapes.
    ``_arr_poly_`` is required
    """

    if ids is None:
        ids = np.arange(len(a)).tolist()
    polys = []
    for i in a:
        p = _arr_poly_(i, sr)  # ---- use _arr_poly
        polys.append(p)
    out = list(zip(polys, ids))
    name = gdb + "\\" + fname
    wkspace = arcpy.env.workspace = 'in_memory'
    arcpy.management.CreateFeatureclass(wkspace, fname, p_type,
                                        spatial_reference=sr)
    arcpy.management.AddField(fname, 'ID_arr', 'LONG')
    with arcpy.da.InsertCursor(fname, ['SHAPE@','ID_arr']) as cur:
        for row in out:
            cur.insertRow(row)
    out_fname = fname + "_mp"
    arcpy.management.Dissolve(fname, out_fname, "ID_arr",
                              multi_part="MULTI_PART",
                              unsplit_lines="DISSOLVE_LINES")
    arcpy.management.CopyFeatures(out_fname, name)
    del cur
    return

Some Results

So just to compare the geometries, I will compare the areas calculated using arcpy and those calculated using the function below.

polys = fc_shapes(in_fc)

areas = [p.area for p in polys]

areas1 = poly_area(a, ids, fr_to)

areas   #   [100.0, 78.0, 155.0, 52.0, 36.0]

areas1  #   array([100.,  78., 155.,  52.,  36.])

Looks good.  The helper function employs my favorite numpy function einsum, to implement the shoelace formula.  A tad overkill for these teeny polygons, but it works blazingly fast for huge point arrays representing real world polygons.

def poly_area(a, ids, fr_to=None):
    """Calculate of a 2D array sliced into sections using an indices of the
    bounds.  ``a`` is created from ``poly_pnts``
    """

    def _e_area(a):
        """mini e_area with a twist, shoelace formula using einsum"""
        x0, y1 = (a.T)[:, 1:]
        x1, y0 = (a.T)[:, :-1]
        e0 = np.einsum('...i,...i->...i', x0, y0)
        e1 = np.einsum('...i,...i->...i', x1, y1)
        return np.nansum((e0-e1)*0.5)
    # ----
    subs = [_e_area(a[f:t]) for f,t in fr_to]
    totals = np.bincount(ids, weights=subs)
    return totals

I could calculate the area and length properties as I construct the arrays, like cursors do so that the property is readily available.  I find if I wanted those properties I would calculate them, and I don't need them floating around unused.

 

poly_area is only one of my helper functions for calculating poly properties.  These will be assembled at some stage and posted on GitHub and the Code Sharing site.

 

 

So why the whole exercise of converting the polygons to arrays in the first place.

Remember, the output is a 2D array of points.  Shifting, rotating, scaling, thinning, anything-ing is done on the whole array at once, not one point at a time.

Also note, that the Polygon and Polyline classes have properties and methods like intersect, union etc.  Simple functions entailing geometry are not shown/available directly from within those classes.  How do you shift a polygon by a finite amount using arcpy? Your homework.

 

Coming soon

The pesky attributes

Part 5 ... The attributes attached to the geometry

 

---- As before, the inputs ----

The polygons that I will be using are shown to the right.

  1. A square, 5 points, first and last duplicates
  2. A lake on an island in a lake...
  3. A multipart with a two different shaped donut holes
  4. The letter 'C'
  5. The letter 'V'

Each part is labelled at the labelpoint rather than the centroid, hence each part gets labelled.

 

 

 

 

 

---- An Alternate Geometry Reconstructor ----
---- Arrays to Poly* Features ----

 

Never use <null> in a table.  To many posts on the forum on how to trap them, find them, replace them.

Always put in a value to represent all conditions.  Too many people use None <null> as the catchall category.  In reality all observations need to be classified exactly, there really is no such thing as <null>.  You either made and observation or you didn't.  If you didn't, your classification scheme should provide a key indicating that.

 

If an observation was made but the phenomenon/parameter/whatever was actually not there, there should be a key for that.  Similar, for 'I forgot', 'Wasn't my job' or whatever other excuses exist.  

 

  • For observations recorded as floating point numbers, that truly yielded 'nothing/None/nadda/zilch', I use 'nan' (np.nan) since it is a recognized number.
  • For text/string observations, I use None since None the object translates to 'None' the text easily in most tables.
  • For time, there is now 'not a time' (NaT), but I don't work with time, preferring to use the string incarnations of those observations
  • For integers... sadly there is no 'Nint'.  You have to provide an actual integer value to represent no value observed, although you desperately tried.   You can use the old school -999, or even 2**8, 2**16 etcetera.
  • For all of the above, anything that is truly doesn't represent an observation where no value was observed, you will have to provide alternatives

 

Making none/null/real nothingness

 

You can add to, or remove from, the list below.  These are some that I use.  I will draw your attention to the NumPy incarnations for integers.  Equivalent values exist for floats.  Any value that ensures that you will take a second look if a calculation looks weird is good.  However if your table contains <nulls> even after my lecture above, this will help mitigate your... stupidity is such a harsh word... but you get my drift

def _make_nulls_(in_fc, int_null=-999):
    """Return null values for a list of fields objects, excluding objectid
    and geometry related fields.  Throw in whatever else you want.

    Parameters
    ----------
    in_flds : list of field objects
        arcpy field objects, use arcpy.ListFields to get a list of featureclass
        fields.
    int_null : integer
        A default to use for integer nulls since there is no ``nan`` equivalent
        Other options include

    >>> np.iinfo(np.int32).min # -2147483648
    >>> np.iinfo(np.int16).min # -32768
    >>> np.iinfo(np.int8).min  # -128

    [i for i in cur.__iter__()]
    [[j if j is not None else -999 for j in i ] for i in cur.__iter__() ]
    """
    nulls = {'Double': np.nan, 'Single': np.nan, 'Float': np.nan,
             'Short': int_null, 'SmallInteger': int_null, 'Long': int_null,
             'Integer': int_null, 'String':str(None), 'Text':str(None),
             'Date': np.datetime64('NaT')}
    #
    desc = arcpy.da.Describe(in_fc)
    if desc['dataType'] != 'FeatureClass':
        print("Only Featureclasses are supported")
        return None, None
    in_flds = desc['fields']
    shp = desc['shapeFieldName']
    good = [f for f in in_flds if f.editable and f.name != shp]
    fld_dict = {f.name: f.type for f in good}
    fld_names = list(fld_dict.keys())
    null_dict = {f: nulls[fld_dict[f]] for f in fld_names}
    # ---- insert the OBJECTID field
    return null_dict, fld_names

My favorite way of getting just the attributes

 

Such nice functions, FeatureClassToNumPyArray, TableToNumPyArray, and back the other way. 

I am sure many of you have explored where it all comes from only to find it all buried in a *.pyd file

import arcpy.da as apd
apd.__file__ # ---- 'C:\\...install path...\\Resources\\ArcPy\\arcpy\\da.py'

# ---- which is actually just imports arcgisscripting
# ---- so import it directly

import arcgisscripting as ags

ags.__file__
# ---- 'C:\\...install path...\\bin\\Python\\envs\\arcgispro-py3\\lib\\
# site-packages\\arcgisscripting.pyd'

 

No bother, since you can pull out data for the attributes nicely accounting for <null> records.

def fc_data(in_fc, verbose=False):
    """Pull all editable attributes from a featureclass tables.  During the
    process, <null> values are changed to an appropriate type.

    Parameters
    ----------
    in_fc : text
        Path to the input featureclass
    verbose : boolean
        Requires ``'prn_rec' in locals().keys()`` in order to set to ``True``.
        ``prn_rec`` is imported from arraytools.frmts
    """

    flds = ['OID@', 'SHAPE@X', 'SHAPE@Y']
    null_dict, fld_names = _make_nulls_(in_fc, int_null=-999)
    fld_names = flds + fld_names
    new_names = ['OID_orig', 'X_centroid', 'Y_centroid']
    a = arcpy.da.FeatureClassToNumPyArray(in_fc, fld_names, skip_nulls=False,
                                          null_value=null_dict)
    a.dtype.names = new_names + fld_names[3:]
    if verbose:
        try: prn_rec(a)  # ** prn_rec imported from arraytools.frmts
        except: print(a)
    return np.asarray(a)

The explorers amongst us, may have discovered a few searchcursor shortcuts

fld_names = ['OBJECTID', 'Long_1', 'Short_1', 'Float_1', 'Double_1', 'Text_1']

cur = arcpy.da.SearchCursor(in_fc, fld_names, explode_to_points=False)

z = cur._as_narray()

Traceback (most recent call last):
File "<ipython-input-220-b4e724f0982b>", line 1, in <module>
z = cur._as_narray()

Sadly, the integer fields with <nulls> bring the whole shortcut down.

The searchcursor actually has enough information in it to create a structured/recarray. 

If you have a clean table with no nulls, the actual calls to _dtype and fields show that you can clearly link cursors and NumPy arrays.  Too bad, the whole integer fix isn't incorporated, but _as_narray and FeatureClassToNumPyArray and TableToNumPyArray yield the same results on a 'clean' dataset.

dir(cur)

[...snip... '_as_narray', '_dtype', 'fields', 'next', 'reset']

 

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

Coming soon

 

  • Work with the geometry... and/or … work with the attributes, then put it all back together.

 

Geometry in NumPy... # 1 

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

Geometry ... Deconstructing poly* features  # 3 

Geometry ... Reconstructing Poly Features # 4