Skip navigation
All People > bixb0012 > Tilting at Globes > 2017 > June
2017

My previous blog post, A Case of Missing Prefixes: ArcGIS ...Geometries, was unplanned.  I started writing this blog post and the next one, and I realized the background information I wanted to include didn't quite fit in either, so I wrote the background information up separately.

 

Over the years, I have found the ArcPy Geometry constructors to be quirky, especially with edge or niche cases like empty geometries.  As as engaged user, who believes providing feedback to developers is important, I have submitted several bugs and enhancement requests over time relating to geometry constructors.  Right off the top, I can't say how many have been implemented, possibly just one, but it wasn't one high up my wish list.

 

A couple years back I submitted the Support Empty Geometries with ArcPy Geometry Class Constructors idea.  Seeing how few users work with empty geometries, I can't say I am surprised with the whopping 2 votes it has received so far.  As David Wynne from Esri points out in the comments, empty geometries are actually supported with the ArcPy Geometry constructors, it is just the syntax is undocumented.  Given empty geometries are supported by the constructors, I think my original idea should be merged under Pythonic ArcPy Geometries since what I am really after is a more intuitive and idiomatic way to create empty geometries using ArcPy Geometry constructors.  I would merge my two ideas if I could, but I can't myself, and I don't want to just delete the idea.

 

Using the ArcPy Polygon class as an example, let's take a look at a variety of ways one might try to construct an empty geometry using the ArcPy Geometry constructors:

>>> ctor_stmts = (
...     '',
...     'None',
...     'arcpy.Array()',
...     'arcpy.Array(None)',
...     'arcpy.Point()',
...     'arcpy.Array([])',
...     'arcpy.Array([None])',
...     'arcpy.Array(arcpy.Array())',
...     'arcpy.Array(arcpy.Array(None))',
...     'arcpy.Array(arcpy.Array([]))'
... )
...
>>> for i, stmt in enumerate(ctor_stmts):
...     stmt = "arcpy.Polygon({})".format(stmt)
...     try:
...         print("{:0>2}. {:<46}: {}".format(i+1, stmt, eval(stmt).WKT))
...     except Exception as err:
...         print("{:0>2}. {:<46}: {}....".format(i+1, stmt, err.message[:35]))
...        
01. arcpy.Polygon()                               : Object: CreateObject cannot create ....
02. arcpy.Polygon(None)                           : Object: CreateObject cannot create ....
03. arcpy.Polygon(arcpy.Array())                  : Object: CreateObject cannot create ....
04. arcpy.Polygon(arcpy.Array(None))              : MULTIPOLYGON EMPTY
05. arcpy.Polygon(arcpy.Point())                  : MULTIPOLYGON EMPTY
06. arcpy.Polygon(arcpy.Array([]))                : Object: CreateObject cannot create ....
07. arcpy.Polygon(arcpy.Array([None]))            : MULTIPOLYGON EMPTY
08. arcpy.Polygon(arcpy.Array(arcpy.Array()))     : MULTIPOLYGON EMPTY
09. arcpy.Polygon(arcpy.Array(arcpy.Array(None))) : MULTIPOLYGON EMPTY
10. arcpy.Polygon(arcpy.Array(arcpy.Array([])))   : MULTIPOLYGON EMPTY
>>>

 

Before sharing any thoughts or comments, I think it is important to refresh ourselves on exactly what the ArcPy Polygon documentation states:

Syntax

 Polygon  (inputs, {spatial_reference}, {has_z}, {has_m})
ParameterExplanationData Type
inputs

The coordinates used to create the object. The data type can be either Point

or Array objects.

Object

 

A few observations before comments:

  • Examples 01 & 02:
    • The documentation states the inputs need to be Point or Array objects, which neither of these examples provide.
  • Examples 03 & 06:
    • These two examples are the same, e.g., repr(arcpy.Array()) == repr(arcpy.Array([]))Since Esri doesn't implement equivalence with ArcPy Arrays, I used the objects' __repr__ to show equality in this case.
    • Calling arcpy.Array() creates an empty ArcPy Array:  <Array []>.  Empty ArcPy Arrays are valid objects, and the documentation states ArcPy Array objects can be passed to the ArcPy Polygon constructor.
  • Examples 04 & 07:
  • Example 05:
    • The documentation states a Point object can be used as an input.
    • Calling arcpy.Point() doesn't create an empty point, like calling arcpy.Array(), instead it creates a 0, 0 point: <Point (0.0, 0.0, #, #)>.  Although not demonstrated above, further testing shows any point works, e.g., arcpy.Polygon(arcpy.Point(10000.5, -4500.3)) also creates an empty multi-polygon.
  • Examples 08 & 10:
    • Similar to Examples 03 & 06 and Examples 04 & 07, these two examples are the same.
    • Calling arcpy.Array(arcpy.Array()) creates an empty ArcPy Array nested within an ArcPy Array.  The ArcPy Array documentation does state an ArcPy Array is a valid input to the ArcPy Array constructor.

 

There are several oddities with the results from the polygon constructor test above.  I think Examples 03 & 06 are most likely to trip users up.  The documentation clearly states ArcPy Array objects are valid inputs, and empty ArcPy Arrays are valid objects, so why does passing an empty ArcPy Array to the ArcPy Polygon constructor generate an error?  Although Examples 04 & 07 work, it does seem odd that an ArcPy Array containing a Python None is somehow substantively different than an empty ArcPy Array when building an empty geometry.

 

Example 05 is a surprise and makes the least sense to me.  The documentation clearly states a Point object is a valid input to the ArcPy Polygon constructor, but what does a polygon with a single point look like?  An empty multi-polygon, according to Esri.  What takes the strangeness up another level for me is that there is no empty ArcPy Point, and that any ArcPy Point that is passed to the ArcPy Polygon constructor creates an empty geometry.

 

Using dis — Disassembler for Python bytecode to look at the bytecode for each of the constructor statements above, it appeared Example 05 would be the highest performing constructor statement.  Running timeit — Measure execution time of small code snippets on the constructor statements showed Example 05 was the quickest at generating empty geometries, not by much, but still the quickest.  This only adds to the mystery of the single-point polygon.

 

Examples 08 & 10 will be touched on in a subsequent blog post, so I won't dive into my thoughts on them here.

 

So how do the results from this test fit into the idea/request/plea for more Pythonic ArcPy Geometries?  For me, the answer is that Example 01 should be the primary way to create ArcPy empty geometries, i.e., passing no arguments to an ArcPy Geometry constructor should generate the equivalent empty geometry type.  For objects that support the concept of emptiness, it is fairly common in Python to have passing no arguments to the constructor create an empty object.  Let's look at some examples:

>>> # Common built-in data structures
>>> list()
[]
>>> set()
set([])
>>> dict()
{}
>>> str()
''
>>>
>>> # Some SciPy data structures now bundled with ArcGIS
>>> import pandas
>>> pandas.DataFrame()
Empty DataFrame
Columns: []
Index: []
>>> import matplotlib.pylab
>>> matplotlib.pylab.plot()
[]
>>>

 

What is more Pythonic for creating an empty geometry, arcpy.Polygon(arcpy.Array(None)) or arcpy.Polygon() ?

 

UPDATE 06/2018:

 

For the sake of completeness, I want to include the results of running the constructor tests on the other ArcPy Geometry classes:

>>> geometry_classes = ["point", "multipoint", "polyline"]
>>> ctor_stmts = (
    '',
    'None',
    'arcpy.Array()',
    'arcpy.Array(None)',
    'arcpy.Point()',
    'arcpy.Array([])',
    'arcpy.Array([None])',
    'arcpy.Array(arcpy.Array())',
    'arcpy.Array(arcpy.Array(None))',
    'arcpy.Array(arcpy.Array([]))'
)
>>> for geometry in geometry_classes:
...     for i, stmt in enumerate(ctor_stmts):
...         stmt = "arcpy.Geometry('{}',{})".format(geometry, stmt)
...         try:
...             print("{:0>2}. {:<60}: {}".format(i+1, stmt, eval(stmt).WKT))
...         except Exception as err:
...             print("{:0>2}. {:<60}: {}....".format(i+1, stmt, err.message[:20]))
...     print("")
...    
01. arcpy.Geometry('point',)                                    : Object: CreateObject....
02. arcpy.Geometry('point',None)                                : Object: CreateObject....
03. arcpy.Geometry('point',arcpy.Array())                       : Object: CreateObject....
04. arcpy.Geometry('point',arcpy.Array(None))                   : Object: CreateObject....
05. arcpy.Geometry('point',arcpy.Point())                       : POINT (0 0)
06. arcpy.Geometry('point',arcpy.Array([]))                     : Object: CreateObject....
07. arcpy.Geometry('point',arcpy.Array([None]))                 : Object: CreateObject....
08. arcpy.Geometry('point',arcpy.Array(arcpy.Array()))          : Object: CreateObject....
09. arcpy.Geometry('point',arcpy.Array(arcpy.Array(None)))      : Object: CreateObject....
10. arcpy.Geometry('point',arcpy.Array(arcpy.Array([])))        : Object: CreateObject....

01. arcpy.Geometry('multipoint',)                               : Object: CreateObject....
02. arcpy.Geometry('multipoint',None)                           : Object: CreateObject....
03. arcpy.Geometry('multipoint',arcpy.Array())                  : Object: CreateObject....
04. arcpy.Geometry('multipoint',arcpy.Array(None))              : MULTIPOINT EMPTY
05. arcpy.Geometry('multipoint',arcpy.Point())                  : MULTIPOINT ((0 0))
06. arcpy.Geometry('multipoint',arcpy.Array([]))                : Object: CreateObject....
07. arcpy.Geometry('multipoint',arcpy.Array([None]))            : MULTIPOINT EMPTY
08. arcpy.Geometry('multipoint',arcpy.Array(arcpy.Array()))     : MULTIPOINT EMPTY
09. arcpy.Geometry('multipoint',arcpy.Array(arcpy.Array(None))) : MULTIPOINT EMPTY
10. arcpy.Geometry('multipoint',arcpy.Array(arcpy.Array([])))   : MULTIPOINT EMPTY

01. arcpy.Geometry('polyline',)                                 : Object: CreateObject....
02. arcpy.Geometry('polyline',None)                             : Object: CreateObject....
03. arcpy.Geometry('polyline',arcpy.Array())                    : Object: CreateObject....
04. arcpy.Geometry('polyline',arcpy.Array(None))                : MULTILINESTRING EMPTY
05. arcpy.Geometry('polyline',arcpy.Point())                    : MULTILINESTRING EMPTY
06. arcpy.Geometry('polyline',arcpy.Array([]))                  : Object: CreateObject....
07. arcpy.Geometry('polyline',arcpy.Array([None]))              : MULTILINESTRING EMPTY
08. arcpy.Geometry('polyline',arcpy.Array(arcpy.Array()))       : MULTILINESTRING EMPTY
09. arcpy.Geometry('polyline',arcpy.Array(arcpy.Array(None)))   : MULTILINESTRING EMPTY
10. arcpy.Geometry('polyline',arcpy.Array(arcpy.Array([])))     : MULTILINESTRING EMPTY

>>>

 

arcpy.Polyline() behaves the same as arcpy.Polygon() while arcpy.Multipoint() has one difference with passing arcpy.Point()

In Semantic Overloading, I touch on one of the main motivations for writing the Tilting at Globes blog:

I believe semantics are important in all aspects of life.  Whether in law, medicine, science, business, information technology, or any other field; having a common language doesn't do much good if there isn't a common understanding of the words that make up the language.  Since languages evolve, maintaining a common understanding of words over time is a continual challenge.

 

Ideally, a common understanding would mean a word or term has a single meaning, and that singular meaning is known and understood by all that use the word or term.  Unfortunately, reality isn't always ideal, and sometimes the context of the word or term plays a large part in its meaning.  For example, how the term "large-scale" is applied to maps can seem counter intuitive to how it is applied to actions, events, and typical objects.  When one understands that map scale applies to the representation of data in the map, i.e., the ratio of distance in the map to distance on the ground, it helps explain why "large-scale" maps cover small geographic areas instead of large geographic areas.  Sure, a large-scale map could cover a large geographic area if printed on an enormous medium, but I am assuming typical print sizes.

 

Moving from maps to geometries, I bring up a case of spatial operators in What's Within:  When (Esri != Clementini) = ? .  In the blog post, I point out how the answer to whether one geometry is "within" another geometry can depend on whom you ask, or more accurately, whose definition of "within" is used to answer the question.  Often times, the qualifiers are either incomplete or left off of the answers entirely, and we end up relying on context to fill in the gaps.

 

In addition to how geometries relate to each other, there are even different contexts for understanding the structure of geometries or geometry types.  Esri's current-generation ArcGIS Geometry Object Model has been around since the release of ArcMap 8.0.  Although the model has evolved some over time, it is not significantly different than when it was released in 1999.  I am not clear on the history of Esri's ST_Geometry storage type, but I know it has been around since at least the ArcGIS 9.x days.  Even operating within a context of Esri geometry types, dropping qualifiers can lead to ambiguity at times.  It can be argued that Esri's REST API represents a 3rd Esri geometry model, but adding that geometry model into this blog post doesn't change the observations and conclusions.

 

Instead of diving into object model diagrams and reference documentation, let's look at an example of how ArcGIS and ST_Geometry types differ in terms of structure.   Below is an image of, and the code to create, four geometries in a SQLite database using Esri's ST_Geometry type.  Each geometry is attributed with both its ST_Geometry and ArcGIS geometry type.

 

ArcMap layout with ST_Geometry sample geometries

 

import arcpy
import os
import sqlite3
from itertools import chain

# Define geometry samples
PCSCode = 3857
geoms = [
    ("poly", "polygon", "POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))"),
    ("multipoly", "multipolygon", ("MULTIPOLYGON(((15 0, 25 0, 25 10, 15 10, 15 0)),"
                                           "((18 13, 18 18, 24 18, 24 13, 18 13)))")),
    ("line", "linestring", "LINESTRING(3 13, 3 18)"),
    ("multiline", "multilinestring", "MULTILINESTRING((6 13, 9 13),(9 14, 9 17))")
]

# Define SQLite database initialization parameters
sqlitedb = r"D:\tmp\test.sqlite" # path to SQLite DB
st_geometry_dll = r"DatabaseSupport\SQLite\Windows32\stgeometry_sqlite.dll"
st_geometry_dll = os.path.join(arcpy.GetInstallInfo()["InstallDir"],
                               st_geometry_dll)
sql_stmts = [[
        "SELECT load_extension('%s','SDE_SQL_funcs_init');" % st_geometry_dll
    ],[
        ("CREATE TABLE %s "
           "(id integer primary key autoincrement not null, "
             "stgeom_type text, arcgis_type text);" % table)
        for table, geom_type, WKT in geoms
    ],[
        ("SELECT AddGeometryColumn "
           "(null, '%s', 'geom', %s, '%s', 2, 'null');" % (table, PCSCode, geom_type))
        for table, geom_type, WKT in geoms
    ],[
        ("INSERT INTO %s "
           "(geom) VALUES (st_geometry('%s', %s));" % (table, WKT, PCSCode))
        for table, geom_type, WKT in geoms
    ],[
        ("UPDATE %s "
           "SET stgeom_type = st_geometrytype(geom);" % table)
        for table, geom_type, WKT in geoms
    ]
]

# Create SQLite database and setup tables
arcpy.CreateSQLiteDatabase_management(sqlitedb, "ST_GEOMETRY")
conn = sqlite3.Connection(sqlitedb)
conn.enable_load_extension(True)
cur = conn.cursor()
for sql in chain(*sql_stmts):
    cur.execute(sql)
conn.commit()
del cur, conn

# Populate ArcGIS geometry type
for table, geom_type, WKT in geoms:
    with arcpy.da.UpdateCursor(
        os.path.join(sqlitedb, table),
        ["shape@", "arcgis_type"]
    ) as cur:
        for shape, _ in cur:
            cur.updateRow([shape, shape.type])
del cur

 

A few observations:

  • ArcGIS has a single polyline geometry type, ST_Geometry has ST_LINESTRING and ST_MULTILINESTRING types.
  • ArcGIS has a single polygon geometry type, ST_Geometry has ST_POLYGON and ST_MULTIPOLYGON types.
  • Whereas ArcGIS has point and multipoint there is no multi- prefix when working with polyline and polygon.

 

If the ST_Geometry types look familiar to those who work with open standards, it isn't coincidence.  The documentation on How is ST_Geometry implemented?—Help | ArcGIS Desktop states in several places that ST_Geometry "is a high-performance storage type that includes ISO- and OGC-compliant SQL access to spatial data."  There are several ISO- and OGC- standards when it comes to geometries, and the "SQL access to spatial data" part of the statement is actually quite important as a qualifier.

 

ISO- and OGC-compliance is fairly broad across various geometry models or storage types, e.g., Microsoft's STGeometry, Oracle's SDO_GEOM,  PostGIS/PostgreSQL ST_Geometry, and others.  The OGC/OpenGIS geometry object model is described in OpenGIS® Implementation Standard for Geographic information - Simple feature access - Part 1: Common architecture.  Instead of embedding a geometry class hierarchy diagram or listing all of the geometry types, I will share there are LineString and MultiLineString types as well as a Polygon and MultiPolygon types.

 

All of this leads to a point, and yes I do have a point, that ArcGIS polyline and polygon types are basically multi-types, just with the multi- prefix missing.  Worse yet, most of the ArcGIS documentation abstracts the user even further from the structure of the ArcGIS Geometry Object Model.  Instead of polylines having paths and polygons having rings, which they do, everything just has "parts."  I touch on this problem of "parts" in The Single Multipart:  Polygons With Interior Boundaries:

The problem with the ArcPy Geometry classes, at least one of them, is that Esri replaced multiple, specific geometry components with a single generic one, the "part."  By abstracting geometries and rings with parts, not only did they deviate from geospatial standards and norms, they created the sideshow-worthy multi-part single-part polygon.

 

Whether it is missing prefixes or the ubiquitous part, understanding the nuances of the implementation and documentation of the ArcGIS Geometry Object Model can help one understand how ArcGIS software, including ArcPy, interacts with other geometry types.