Point class in arcpy... why the limits?

DanPatterson_Retired · ‎03-03-2015

Create four points with slightly different properties to test other properties.

>>> pnt0 = arcpy.Point(X=1,Y=2,Z=3,M=4,ID=0) # a point
>>> pnt1 = arcpy.Point(X=1,Y=2,Z=1,M=0,ID=1) # different Z,M and ID
>>> pnt2 = arcpy.Point(X=3,Y=4,Z=2,M=2,ID=2) # different X,Y,Z,M and ID
>>> pnt3 = arcpy.Point(X=1,Y=2,Z=3,M=4,ID=3) # clone pnt0, changing ID only
‍‍‍‍

Perform some equality checks to test what is comparable to pnt0. There is some fancy attribute collection and print formatting.

>>> all_pnts = [pnt0,pnt1,pnt2,pnt3]
>>> equality = [[getattr(pnt0,i) == getattr(p,i) for i in props] for p in all_pnts]
>>> for i in range(len(equality)):
  ... print('pnt0 vs pnt{}: {}'.format(i, equality[i]))
  ...

pnt0 vs pnt0: [True, True, True, True, True]
pnt0 vs pnt1: [True, True, False, False, False]
pnt0 vs pnt2: [False, False, False, False, False]
pnt0 vs pnt3: [True, True, True, True, False]
>>>
‍‍‍‍‍‍‍‍‍‍‍

It is good, that pnt0 is identical to itself, and pnt3 only differs by the ID number. All the other comparisons identify what attributes in the properties list (prop) differ. Why not simplify things then and use the built-in methods? Well...because...

>>> pnt0 == pnt3
False
>>> pnt0.equals(pnt3)
True
>>>
‍‍‍‍‍

The first equality check ( == ) is comparing objects and all its properties, whilst the second

check doesn’t consider the ID value. To confirm this consider the following tests.

I would expect this behaviour

real_pnt = arcpy.Point(0,0,1,1,0)
>>> real_pnt2 = arcpy.Point(0,0,1,1,0)
>>> real_pnt == real_pnt
True
>>> real_pnt == real_pnt2
False
‍‍‍‍‍‍

not this

>>> real_pnt.equals(real_pnt2)
True‍‍

>>> a = arcpy.Point(X=1,Y=2,Z=3,M=4)
>>> b = arcpy.Point(X=1,Y=2,Z=3,M=0)
>>> c = arcpy.Point(X=1,Y=2,Z=0,M=0)
>>> a == b, a == c, b == c
(False, False, False)
>>> a.equals(b), a.equals(c),  b.equals(c)
(True, True, True)
>>>
‍‍‍‍‍‍‍‍

Shades of George Orwell...

"...all Points are equal but some Points are more equal than others..."

Disappointing if you require comparisons on all point properties.

Well PointGeometry is no better, since it depends upon the points that are used to construct it. Consider the equality check again, using PointGeometry objects instead of Point objects.

>>> pg_a = arcpy.PointGeometry(a)
>>> pg_b = arcpy.PointGeometry(b)
>>> pg_c = arcpy.PointGeometry(c)
>>> 
>>> pg_a == pg_b, pg_a == pg_c, pg_b == pg_c
(False, False, False)
>>> pg_a.equals(pg_b), pg_a.equals(pg_c), pg_b.equals(pg_c)
(True, True, True)
>>> 
‍‍‍‍‍‍‍‍‍

Any comments as to why the Arcpy Point and PointGeometry class are so feable? Is the implementation in ArcObjects better? Or should I go back to using my own geometry classes.

I will be producing a summary on my blog after my musings about the various geometry classes is complete but my impressions are somewhat muted by the seemingly poor implementation or exposure to useful properties and methods.

NeilAyres · ‎03-03-2015

Aren't you being a little unfair Dan?

The

geom1.equals(geom2)

is an arcpy geometry method, and the help under arcpy Point clearly states that this is a 2D comparison only.

Whereas "==" is a python equivalence test. They do different things.

DanPatterson_Retired · ‎03-03-2015

Nope...I know it clearly states it... it doesn't make it right. Also, working with Z and/or M geometry is far simpler in other packages. Also..

>>> import arcpy
>>> null_pnt = arcpy.Point()
>>> null_pnt
<Point (0.0, 0.0, #, #)>
>>>
‍‍‍‍‍

Really? isn't 0 for X and Y supposed to be valid coordinates?

or this

>>> 
>>> real_pnt = arcpy.Point(0,0,1,1,0)
>>> real_pnt
<Point (0.0, 0.0, 1.0, 1.0)>
>>> null_pnt == real_pnt
False
>>> null_pnt.equals(real_pnt)
True
>>>
‍‍‍‍‍‍‍‍‍

Not good

JoshuaBixby · ‎03-03-2015

In terms of null points, I completely agree with you that arcpy.Point() should not return the <Point (0.0, 0.0, #, #)> that it does. I think Shapely/GEOS/JTS handles the situation much better:

>>> from shapely.geometry import Point
>>> null_pnt = Point()
>>> print null_pnt
GEOMETRYCOLLECTION EMPTY

With an empty geometry collection, comparing a null/empty point to a valid regular point with 0, 0 coordinates behaves the way one would suspect:

>>> real_pnt = Point(0, 0, 1)
>>> print real_pnt
POINT Z (0 0 1)
>>> null_pnt.equals(real_pnt)
False

I think an empty geometry collection is the way to go; unfortunately, Esri hasn't implemented a geometry collection type with ArcPy. Off the top, four different ideas come to mind how to handle the situation without a geometry collection:

have arcpy.Point() throw an error saying the coordinates aren't valid
implement an empty point concept and return an empty point
return a None object
return a valid point using the 0, 0 coordinates

Out of the four options above, I think #4 (the existing condition) is the least desirable from a user/developer perspective.

JoshuaBixby · ‎03-03-2015

I agree that Python object equivalence/equality likely isn't the best approach for comparing spatial objects primarily because I haven't seen any geospatial Python packages where the various geometry classes' __eq__() methods are spatially aware. With ArcPy, they are not.

>>> pl1 = arcpy.Polyline(
...     arcpy.Array(
...         [arcpy.Point(0, 0), arcpy.Point(1, 1)]
...     )
... )
... 
>>> pl2 = arcpy.Polyline(
...     arcpy.Array(
...         [arcpy.Point(0, 0), arcpy.Point(0.5, 0.5), arcpy.Point(1, 1)]
...     )
... )
... 
>>> pl1 == pl2
False
>>> pl1.equals(pl2)
True

For points, Python object equality works in some situations because of how simple points are spatially. Once we move into polylines and polygons, and possibly multipoints, using Python object equality falls short.

JoshuaBixby · ‎03-03-2015

How ArcPy is handling methods for testing spatial relations between geometric objects, e.g., equals, is consistent with the OGC Simple Feature Access - Part 1: Common Architecture use of Z and M coordinate values:

6.1.2.5 Use of Z and M coordinate values
A Point value may include a z coordinate value....
A Point value may include an m coordinate value....
Observer methods returning Point values include z and m coordinate values when they are present.
Spatial operations work in the "map geometry" of the data and will therefore not reflect z or m values in calculations (e.g., Equals, Length) or in generation of new geometry values (e.g., Buffer, ConvexHull, Intersection). This is done by projecting the geometric objects onto the horizontal plane to obtain a "footprint" or "shadow" of the objects for the purposed of map calculations. In other words, it is possible to store and obtain z (and m) coordinate values but they are ignored in all other operations which are based on map geometries. Implementations are free to include true 3D geometric operations, but should be consistent with ISO 19107.

The operative statement for this discussion is really the last one, i.e., implementations are free to include true 3D [or 4D] geometric operations, but I will circle back to this in a minute.

Other products that implement the OGC Simple Feature Access standards behave the same way as ArcPy in this instance. Looking at Shapely, which is based on GEOS and JTS libraries, and grabbing a few example points from the original code snippets:

>>> from shapely.geometry import Point
>>> a = Point(1, 2, 3)
>>> print a
POINT Z (1 2 3)
>>> c = Point(1, 2, 0)
>>> print c
POINT Z (1 2 0)
>>> a == c
False
>>> a.equals(c)
True

Even looking at how MS SQL Server implements its Geometry data type:

>>>Declare @a geometry, @c geometry
>>>SET @a = 'POINT(1 2 3 4)'
>>>PRINT @a.ToString()
POINT (1 2 3 4)
>>>SET @c = 'POINT(1 2 0 0)'
>>>PRINT @c.ToString()
POINT (1 2 0 0)
>>>PRINT @a.STEquals(@c)
1

I think the behavior being seen in ArcPy is typical, some might even argue expected given the Simple Feature specifications/standards. That said, Esri could implement OGC Feature Geometry/ ISO 19107 Geographic information -- Spatial schema in ArcPy, but even that would not necessarily address all of the issues here because Esri has created a 5D point model (x, Y, Z, M, ID).

Has Esri failed or dropped the ball on this one? It could be argued, but their existing implementation could also be defended on some levels. The Simple Features specification sets the bar pretty low, but Esri also hasn't tried to raise it for themselves. There is nothing in the specifications that prevents them from adding more robust functionality. Think of how Microsoft handled the situation, i.e., they have OGC methods and Extended methods so they can remain compliant but add extra functionality.

DanPatterson_Retired · ‎03-03-2015

This is the problem I am finding, since I now work quite a bit with numpy arrays for a variety of reasons and handling 3D and above coordinates is a breeze the problem getting things out to and back from numpy. As an example, I tried to create a cube which I can do quite readily in a variety of packages that I use. Try to create a polygon that is perfectly vertically aligned...not extruded...but oriented 90 degrees to the axis. No matter what output format I tried, I got nothing and I don't want to use TINs etc.

I even looked at GeoJson etc and Arc kills the m coordinate even though the files are z and m aware, if you have X,Y and Z, I need the M to use homogenous coordinates to do my stuff.. Oh well, any more thoughts or workarounds are appreciated.

JoshuaBixby · ‎03-04-2015

Overall, I think the ArcPy support for Z and M coordinates is poor, quite poor. I wish I had a better response, but I just hit walls/bugs/design limitations around every corner.

I quickly revisited arcpy.AsShape and arcpy.FromWKT. Although the Help for FromWKT doesn't state anything one way or another about Z or M coordinates, the function throws an error when you try to use Z or M coordinates. I have been able to use arcpy.ArcSDESQLExecute to pass Z and M WKT coordinates to SQL Server Express, basically bypass ArcPy and relying on the SQL Server engine to handle creating the shapes correctly. Even though it has worked for me, and I can view the results in ArcMap, all of the regular and Data Access cursors seem to do nothing with them. It might just be ArcPy cursors can't cope with Z or M coordinates, need to investigate more.

Good luck, let us know if you have any success.

UPDATE: I had to strike through some of my earlier comments. On the WKT comment, it is still the case that the Help doesn't explicitly mention Z or M coordinates; however, arcpy.FromWKT was throwing an error because I wasn't giving it OGC compliant WKT representations. I am used to MS SQL Server's WKT implementation, which is slightly different than the OGC specifications, so what I would pass SQL Server was failing with ArcPy.

On the ArcPy cursors not working with ZM polygons created through SQL Server and arcpy.ArcSDESQLExecute, I had a goofy projection discrepancy that was causing the problems. Once I went back and straightened out the projections, I was able to use cursors with the newly created polygons.

DanPatterson_Retired · ‎03-04-2015

thanks Joshua... I am currently sub-typing ndarrays and recarrays to make what I do more user friendly when working with geometry objects. I will have to worry about visualizing the outputs in ArcMap in true 3D and when it uses Z values to discern two poly* features that share common boundaries but differ in Z (either for multipart formation without shape dissolve, or for rendering in 3D). There has to be something in arcobjects. I can't believe that all they have is extrusion and base heights as an option. What about City Engine? is anyone using that?

I will post, when I have more. I hope anyone with further info will chime in.

JoshuaBixby · ‎03-05-2015

Since ArcGIS Pro is an integrated 2D/3D application, it will be interesting to see if Esri extends the existing ArcPy Geometry classes or introduces new ones that support 3D differently.