What's Within: When Esri + Clementini = User Beware

[JoshuaBixby]

This is the second in a two-part series on the risks to software users of poor documentation; specifically, the confusion and unexpected results that come from weakly documented spatial operators in GIS software.  The first part in the series looks at how inconsistent and incomplete documentation requires users to guess how spatial operators are implemented.  The second part in the series looks at the inconsistent results that arise from mixing different implementations of spatial operators.

One of the many empty geometry bugs I have submitted recently got put on the front burner this week when I noticed Esri updated the status to "Closed:  Will Not be Addressed."  What has ensued is a discussion that is still unfolding over expected behavior versus unexpected results.  Coincidentally, the issue can be framed neatly within the discussion and examples from the first post in this series (What's Within:  When (Esri != Clementini) = ?).

The first post in this series discussed how there is no singular definition of spatial relations for geometry libraries and geospatial applications, and how the Dimensionally Extended 9 Intersection Model (DE-9IM) became the prevailing 2D definition after inclusion in the OpenGIS Implementation Specification for Geographic information - Simple feature access - Part 1: Co....  The post focused on how incomplete documentation of spatial operators forces users to guess, and likely incorrectly at times, how software works instead of knowing how the software should work.  As unfortunate for users as incomplete documentation can be, mixing different definitions of spatial relations within the same spatial operator is much worse, and that appears to be what Esri has done.

Borrowing from the examples in the previous post, let's look at three features and the ArcPy Geometry.within() method.

>>> #Create square polygon
>>> polygon = arcpy.FromWKT('POLYGON((0 0, 3 0, 3 3, 0 3, 0 0))')
>>> #Create line that completely lies on polygon boundary
>>> line = arcpy.FromWKT('LINESTRING(1 0, 2 0)')
>>> #Create empty line
>>> line_empty = arcpy.FromWKT('LINESTRING EMPTY')
>>>
>>> #Test whether line is within polygon
>>> line.within(polygon)
False
>>> #Test whether line_empty is within polygon
>>> line_empty.within(polygon)
True

For comparative purposes, let's run the same two tests from Lines 9 and 12 using Esri’s ST_Geometry in Oracle:

SQL> --Test whether line is within polygon
SQL> SELECT sde.st_within(sde.st_geomfromtext('LINESTRING(1 0, 2 0)', 0),
                          sde.st_geomfromtext('POLYGON((0 0, 3 0, 3 3, 0 3, 0 0))', 0))
       FROM dual;
SDE.ST_WITHIN(SDE.ST_GEOMFROMTEXT('LINESTRING(10,20)',0),SDE.ST_GEOMFROMTEXT('PO'
--------------------------------------------------------------------------------
                                                                               0
SQL> --Test whether line_empty is within polygon
SQL> SELECT sde.st_within(sde.st_geomfromtext('LINESTRING EMPTY', 0),
                          sde.st_geomfromtext('POLYGON((0 0, 3 0, 3 3, 0 3, 0 0))', 0))
       FROM dual;
SDE.ST_WITHIN(SDE.ST_GEOMFROMTEXT('LINESTRINGEMPTY',0),SDE.ST_GEOMFROMTEXT('POLY'
--------------------------------------------------------------------------------
                                                                               0

Interesting, the results from using ST_Geometry differ from using ArcPy.  Knowing that ST_Geometry functions are OGC simple feature access and SQL compliant, which means they implement DE-9IM, the results from ST_Geometry in Oracle are expected because a line solely on the boundary of a polygon is not considered within the polygon, and an empty geometry cannot be within another geometry.

When looking at ArcPy results, Line 10 is correct only if the ArcPy Geometry Classes implement Clementini's definition since we saw in the first post of this series that Esri's definition of Within is True in this situation.  Unfortunately, the ArcPy documentation doesn’t state one way or another which definition it is implementing.  The result on Line 13 is correct only if the ArcPy Geometry Classes implement Esri's definition because the Clementini definition does not allow for an empty geometry to be within another geometry.  Clear as mud, right?

Esri's stance, to date, is that everything is working as designed.  What?!  If that is the case, Esri is implicitly acknowledging they are implementing different definitions of a spatial relation within the same spatial operator, it just depends what geometries you pass it!!  Between the closed source code, incomplete documentation, and seemingly arbitrary implementation; ArcPy Geometry Classes are use at your own risk.  Caveat utilitor.