Call for More Pythonic ArcPy Geometries

Discussion created by bixb0012

on Jul 14, 2015
Latest reply on Feb 19, 2016 by bixb0012

Like • Show 0 Likes0
Comment • 17

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.

Outcomes

17 Replies

xander_bakker Jul 14, 2015 3:47 PM
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...
Like • Show 0 Likes0
Actions
- Dan_Patterson @ xander_bakker on Jul 14, 2015 3:50 PM
  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)
  2 people found this helpful
  Like • Show 0 Likes0
  Actions
  - xander_bakker @ Dan_Patterson on Jul 14, 2015 4:03 PM
    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))
    Like • Show 0 Likes0
    Actions
    - Dan_Patterson @ xander_bakker on Jul 14, 2015 4:52 PM
      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.
      1 person found this helpful
      Like • Show 0 Likes0
      Actions
      - bixb0012 @ Dan_Patterson on Jul 15, 2015 10:19 AM
        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.
        Like • Show 3 Likes3
        Actions
        Dan_Patterson @ bixb0012 on Jul 15, 2015 1:33 PM
        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.
        Like • Show 1 Like1
        Actions
    - bixb0012 @ xander_bakker on Jul 16, 2015 3:07 PM
      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.
      Like • Show 0 Likes0
      Actions
Dan_Patterson Jul 16, 2015 3:22 PM
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.
Like • Show 1 Like1
Actions
- bixb0012 @ Dan_Patterson on Jul 18, 2015 7:39 PM
  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.
  3 people found this helpful
  Like • Show 0 Likes0
  Actions
  - Dan_Patterson @ bixb0012 on Jul 19, 2015 4:06 AM
    you should be able to create both. The island version if easily doable but the multipart by attribute is another story. Xander pointed that out to me in a thread (which I can't track down now), that you can make valid multipart polygons if they share no more than a single point but if they share one line then the intervening line gets removed and adjacent polygons get combined EVEN if their Z or M values differ. Picture a series of squares in a row, with increasing Z values but each sharing a common line (or two), there is no way to make that a multipolygon because arcpy and its underlying code only uses X and Y in tests for equality. Also...try to make a vertical polygon for display in 3D...you can't and I have to rely on other software for display. And forget getting distracted and trying to construct a tessaract...Tesseract - Wikipedia, the free encyclopedia
    
    I would have hoped that with the push to 3D that the geometry engine be improved, but I have seen no reference that the new Pro offers much in that...It must lie in City Engine, but I have had no time to check
    1 person found this helpful
    Like • Show 0 Likes0
    Actions
    - xander_bakker @ Dan_Patterson on Jul 19, 2015 6:36 PM
      Hi Dan Patterson,
      
      I Think you are referring to this question I posted earlier...Create a multipart polygon with arcpy
      Interesting thing the Tesseract - Wikipedia, the free encyclopedia , I wonder how far I'll get...
      
      Maybe this will help... tesserpy 1.1.1 : Python Package Index
      1 person found this helpful
      Like • Show 0 Likes0
      Actions
bixb0012 Feb 12, 2016 1:34 PM
Rumor has it that Geometry indexing and slicing is coming at ArcGIS 10.4.1. I am still trying to track down specifics, but what I have heard so far is at least a glimmer of hope.
Like • Show 1 Like1
Actions
- Dan_Patterson @ bixb0012 on Feb 12, 2016 2:09 PM
  keep us posted!!! hopefully before my sabbatical ends
  Like • Show 0 Likes0
  Actions
  - xander_bakker @ Dan_Patterson on Feb 19, 2016 1:32 PM
    This is what's mentioned on the What's new (in ArcMap ) page:
    What's new in ArcMap—ArcGIS Help | ArcGIS for Desktop
    
    Python and ArcPy
    ArcGIS 10.4 has been upgraded to include Python 2.7.10. Additional third-party libraries including SciPy, pandas, Sympy, and nose have been included, and existing third-party libraries including NumPy and matplotlib have been upgraded to more current releases.
    The GetPackageInfo function has been added to return a dictionary about packages.
    
    I'm struggling here at the client site to update the license manager to 10.4, but I have to file a change request and a lot of related stuff... for 5 minutes work... so no way I can test to see if there is more news on the geometries...
    2 people found this helpful
    Like • Show 0 Likes0
    Actions
    - Dan_Patterson @ xander_bakker on Feb 19, 2016 2:16 PM
      Good... some progress... but I can't figure out why they didn't upgrade to at least 3.4.* since ArGIS Pro uses it. So all the SciPy stack modules are for 2.7.10 I presume. I have happily been using the 3.4.* stack for a while, the differences between numpy 1.8 and 1.9 (or 1.10) are subtle, but there have been a few changes that are worth it... I will post more on the diffrence in the Numpy Repository when I get some time.
      Like • Show 0 Likes0
      Actions
      - bixb0012 @ Dan_Patterson on Feb 19, 2016 2:57 PM
        It would be great if Esri gave an indication of when, or if, they will move ArcGIS for Desktop (or is it ArcMap) to Python 3.x. The deeper I get into my transition to 3.x, the more little things pile up that are annoyances going back to 2.7.x.
        Like • Show 0 Likes0
        Actions
    - bixb0012 @ xander_bakker on Feb 19, 2016 3:00 PM
      From what I have heard, the geometry changes are still baking, didn't quite rise in time for the 10.4.0 release.
      
      I hope indexing and slicing bake enough for the 10.4.1 release. I would like to see indexing and slicing return other geometry objects, but I get the impression they will return ArcPy arrays much like getPart.
      Like • Show 0 Likes0
      Actions