Call for More Pythonic ArcPy Geometries

[JoshuaBixby]

I have been working with Django, specifically GeoDjango, on and off for a few weeks, and I must say it has opened my eyes.  When working with the GeoDjango GEOS API, geometries work quite a bit different than in ArcGIS.

Geometries are Pythonic

GEOSGeometry objects are ‘Pythonic’, in other words components may be accessed, modified, and iterated over using standard Python conventions.

Since the idea is already implemented and documented, I won't dive into the specifics here of how true or more Pythonic geometries work.  Given that ArcPy is a Python site package, I think ArcPy geometries should be a bit more robust in terms of using standard Python conventions to interact with them.  This is one area where ArcGIS Pro has been a big let down, i.e.,  ArcGIS Pro was not used as an opportunity to re-invent/re-invigorate the ArcPy site package.

I have posted the idea on ArcGIS Ideas (Pythonic ArcPy Geometries) and opened an Esri Support Enhancement Request (ENH-000089037: Make ArcPy Geometries truly Pythonic) if anyone wants to vote for it or attach their customer number to the idea.

[XanderBakker]

Interesting thought... I just promoted the idea. I wonder what Esri will respond. For the current arcpy (ArcMap) site package I can image that the dependencies to ArcObjects are a huge constraint. One could think that the compatibility of scripts created in 10.x with ArcGIS Pro may have been a consideration, but there will certainly be more reasons for this.

So, Joshua Bixby , what have you been up to with GeoDjango? It is always interesting to hear what other people are doing...

[DanPatterson_Retired]

To add to Xander's requests... some cool examples would be nice...also

Does it answer some of my thread posts reported previously regarding null geometry, outer/inner ring construction (both are clockwise... for example).  And how well are existing long-standing data constructs accessable (ie shapefiles, derivatives in numpy array format etc etc, integration with NumPy and/or SciPy)

[XanderBakker]

As far as I can see from the documentation (GEOS API | Django documentation | Django ) the inner and outer rings are both clockwise:

ext_coords = ((0, 0), (0, 1), (1, 1), (1, 0), (0, 0))
int_coords = ((0.4, 0.4), (0.4, 0.6), (0.6, 0.6), (0.6, 0.4), (0.4, 0.4))
poly = Polygon(ext_coords, int_coords)
poly = Polygon(LinearRing(ext_coords), LinearRing(int_coords))

[DanPatterson_Retired]

​Xander....I saw that but I didn't see the equivalent for a null point

None isn't...nor is 0 or 1 ... more explorations into geometry

where Joshua made reference to Shapely and reference was made to comparable structures in GeoJson, WKT and WKB etc

Then I question how attribute and or geometry deals with the whole concept of 'nothing' (I know his comment link deals with pythonic geometry, but pythonic attributes will be an issues at some time in arcpy)

NumPy Snippets # 6 .... much ado about nothing ... NaN stuff

I shall have to get distracted again, should this be a useful path to pursue since I am currently using a combination of NumPy, my own Point class and arcpy to deal with some of the things that I am doing.  If you get to UC, and find out anything about arcpy plans from 'those in the know' I would be interested in hearing, since ArcGIS Pro offers me little in progress.

[JoshuaBixby]

Xander Bakker​ and Dan Patterson​, great comments and dialogue, thanks for chiming in.  I will elaborate a bit more, but will add comments at the root level to minimize how deep the nested comments go.

Overall, like all software packages/applications/libraries, GeoDjango isn't perfect.  That said, there are some things it does well in terms of implementing the GEOS library and how one interacts with geometry objects.

I think the NumPy integration in ArcGIS/ArcPy is big, a really positive step, but I also see it existing outside of the ArcPy Geometry classes.  Maybe the ArcPy Geometry classes are used under the hood for part of the transition between geospatial data store and NumPy array, but the user doesn't directly interact with ArcPy Geometry classes.

[DanPatterson_Retired]

Which begs the question, why doesn't esri provide more documentation on how to use the association?

In this thread Numpy Snippets # 3 ... Phish_Nyet ... creating sampling grids using numpy and arcpy

I made extensive use of NumPy to produce polygon geometry because I couldn't access what I needed through arcpy and I had no intention of using ArcObjects.  The existing arcpy geometry toolset ... Geometry—ArcPy Classes | ArcGIS for Professionals  ... is quite extensive, however, simple geometric operations like translation, rotation and scaling have long been hidden in a dense access layer since the days of Avenue.  And why in the NumPyArrayToFeatureClass documentation is the only output geometry 'apparently' a point featureclass...tsk.

[JoshuaBixby]

The Polygon constructor is fairly lenient in terms of ring direction/orientation, but the GEOS Polygon object has a normalize() method that will convert the rings to normal/canonical form in place:

>>> ext_coords_cw = ((0, 0), (0, 1), (1, 1), (1, 0), (0, 0))
>>> ext_coords_ccw = ((0, 0), (1, 0), (1, 1), (0, 1), (0, 0)) 
>>> int_coords_cw = ((0.4, 0.4), (0.4, 0.6), (0.6, 0.6), (0.6, 0.4), (0.4, 0.4))
>>> poly_cw = Polygon(ext_coords_cw, int_coords_cw)
>>> poly_ccw = Polygon(ext_coords_ccw, int_coords_cw)
>>> poly_cw.coords
(((0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0), (0.0, 0.0)), ((0.4, 0.4), (0.4, 0.6), (0.6, 0.6), (0.6, 0.4), (0.4, 0.4)))
>>> poly_ccw.coords
(((0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (0.0, 0.0)), ((0.4, 0.4), (0.4, 0.6), (0.6, 0.6), (0.6, 0.4), (0.4, 0.4)))
>>> poly_cw.normalize()
0
>>> poly_cw.coords
(((0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0), (0.0, 0.0)), ((0.4, 0.4), (0.6, 0.4), (0.6, 0.6), (0.4, 0.6), (0.4, 0.4)))
>>> poly_ccw.normalize()
0
>>> poly_ccw.coords
(((0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0), (0.0, 0.0)), ((0.4, 0.4), (0.6, 0.4), (0.6, 0.6), (0.4, 0.6), (0.4, 0.4)))

The GEOS library appears to subscribe to a clockwise outer ring rule, which then makes the inner rings counterclock wise.

[DanPatterson_Retired]

Will have to give it a try on Island in a Lake on an Island in a Lake on an Island : Image of the Day to check the simplicity of the polygon construct, but it is nice that the ring construct is the same and will work nicely with numpy and other libraries.

[JoshuaBixby]

You piqued my interest.  Are you interested in making an invalid Polygon, valid MultiPolygon, or both when you mention wanting to try the island in a lake on an island... example?  The GeoDjango GEOS API will let you make a Polygon with nested holes, but it will also tell you it isn't valid because of those nested holes if you call the valid_reason method.  That is yet another thing lacking from the ArcPy Geometry classes that exists in Shaply and GeoDjango, i.e., the ability to check for valid geometries by simply calling a method or function on the geometry instead of running a geoprocessing tool that has to create an output table with results.