topic Re: Point class in arcpy... why the limits? in Python Questions

Point class in arcpy... why the limits?

DanPatterson_Retired — Sat, 11 Dec 2021 14:54:28 GMT

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.

Re: Point class in arcpy... why the limits?

NeilAyres — Tue, 03 Mar 2015 13:14:39 GMT

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.

Re: Point class in arcpy... why the limits?

DanPatterson_Retired — Sat, 11 Dec 2021 14:54:31 GMT

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

Re: Point class in arcpy... why the limits?

JoshuaBixby — Sat, 11 Dec 2021 14:54:34 GMT

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.

Re: Point class in arcpy... why the limits?

JoshuaBixby — Sat, 11 Dec 2021 14:54:36 GMT

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.

Re: Point class in arcpy... why the limits?

JoshuaBixby — Sat, 11 Dec 2021 14:54:39 GMT

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.

Re: Point class in arcpy... why the limits?

DanPatterson_Retired — Wed, 04 Mar 2015 00:33:35 GMT

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.

Re: Point class in arcpy... why the limits?

JoshuaBixby — Wed, 04 Mar 2015 20:07:42 GMT

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.

Re: Point class in arcpy... why the limits?

DanPatterson_Retired — Wed, 04 Mar 2015 22:44:40 GMT

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.

Re: Point class in arcpy... why the limits?

JoshuaBixby — Thu, 05 Mar 2015 17:22:58 GMT

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.

Re: Point class in arcpy... why the limits?

ElleryChan — Thu, 05 Mar 2015 20:00:50 GMT

It seems important to keep in mind that unless arcpy.Point defines a custom __eq__ or __cmp__ operator, then comparisons like pointA == pointB will only return True if pointA and pointB are references to the same object. Apparently the arcpy developers decided to stick with that behavior and define their own equals() method for geospatial comparisons.

However, IMHO your questioning of the whether the semantics of Point.equals() are adequate is a valid question.

Re: Point class in arcpy... why the limits?

DanPatterson_Retired — Sun, 10 May 2015 00:42:12 GMT

from Joshua's post....

Since ArcGIS Pro is an 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.

any news from the ArcGIS Pro users? or is its incarnation not added to the geometry capabilities much?

Re: Point class in arcpy... why the limits?

curtvprice — Sun, 10 May 2015 04:08:54 GMT

> any news from the ArcGIS Pro users? or is its incarnation not added to the geometry capabilities much?

The help for pro looks awfully similar.

Point—ArcPy Classes | ArcGIS for Professionals

PointGeometry—ArcPy Classes | ArcGIS for Professionals

Re: Point class in arcpy... why the limits?

JoshuaBixby — Mon, 11 May 2015 02:39:09 GMT

Beyond the ArcPy Mapping module (arcpy.mp), which had to get updated quite a bit because of how different the ArcGIS Pro UI is compared to ArcGIS Desktop, it is disappointing all of the other modules and base classes have been updated and cleaned up so little.

Esri took the release of ArcGIS Pro, a brand new application, to migrate ArcPy to Python 3.x, but I would have liked to see the migration to Python 3.x include cleaning up some of the detritus. For example, why keep the iterator.next() methods around with the cursors in ArcGIS Pro, especially the arcpy.da cursors? PEP 3114 was accepted back in 2007 as an acknowledgement that calling for the next item of an iterator in Python 2.x was/is not consistent with how other magic/dunder methods are handled in the language.

Esri did add dunder next methods in ArcGIS Pro because they had to or none of the cursors would work, but they didn't drop the explicit next method. Some might say including the explicit next method allows for a smoother transition from ArcGIS Desktop and Python 2.x to ArcGIS Pro and Python 3.x, but that argument is rather facile. The major overhaul to ArcPy Mapping already broke the smooth transition argument, and keeping an explicit next method with the cursors just promotes poor Python programming, especially for those unfortunate souls who decide to learn Python through ArcPy.

Re: Point class in arcpy... why the limits?

DanPatterson_Retired — Mon, 11 May 2015 06:53:19 GMT

I thank all for the comments. I have 'Helpful' -ed those that provided some context to the questions. In my estimation here are some conclusions:

arcpy may be the holdup and this would include cases such as:
- a spatial reference object is required 'precise' geometry to be created within the confines general possibilities (many examples of errors in outputs because a spatial reference not used...workaround...use one). Highest precision geometry should be the default. I erroneously assumed that it would be a fundamental tenant of a GIS platform.
- equality checks in geometry are 'flawed' at worst or 'incomplete' to be supported. For example, although Z values are included in geometry they are not used in shape geometry checks to see if two geometries differ ... ie (1,1,1 == 1,1,0 in terms of X,Y,Z values... I did not bother to go further to check whether X,Y,M exhibited the same behaviour)
- since equality checks can be called into question ... I presume other Boolean checks, might exhibit similar behaviour... ie. limited to 2D
- dissolving two adjacent shapes based upon an attribute removes the dividing geometry...this should be an option not the default. What if I wanted to show adjacent equal land classes yet retain the boundary as in the case of disjoint multipart shapes.
Given the above, effort seems to have been placed on presentation output and dealing with handling of the mapping module in advance to those dealing with analysis .... disappointing to me
The enhanced capabilities and integration of 2D and 3D into one 'home' seems moot. I want to be able to form a cube or tesseract given coordinates with respect to a reference plane. We appear to be stuck with extrusion to fill in the sides of a cube between two suspended planes since equality checks
digitizing either onscreen or on paper has not evolved, in my estimation....perhaps clouded by students complaining ...'so much easier in CAD...can we use it for the exercise? .... But in the end you will be able to get your frustratingly produced map out in a tiled pdf format using the mapping module

I applaud efforts being made and that the organization is being drawn in different directions and that the biggest return on investment will dictate what 'barking dog' gets fed first .... but come on ... what is wrong with marketing Arc* for its geometric capabilities and quality. What is also wrong with making it a showcase feature.

Some things on my wish list and for my sabbatical:

a Selection tool that reports feature geometry properties on-the-fly
an Add Data tool that does more than report an undefined coordinate system. It should report that the coordinate systems are mixed and give you the option of physically projecting to one coordinate system ... it would sure get rid of a lot of questions on GeoNet
disable any tool that currently enables you to carry out analysis using Geographic Coordinates ... pretty well the whole Spatial Analyst extension in my estimation
I will leave the trailing ellipse for others to ponder
...

Re: Point class in arcpy... why the limits?

JoshuaBixby — Mon, 11 May 2015 15:14:39 GMT

Couldn't agree more. A robust geometry library and implementation should be seen as the foundation of a Geographic Information System, not as some bin item or nice-to-have.

Esri's lack of and inconsistency with supporting empty geometries is another case in point. I opened five or more bugs in the past two weeks on this topic, and I am not battling Esri to justify cleaning up the inconsistencies and adding better support. The Romans did just fine without a zero, right? Now where is my abacus.

Re: Point class in arcpy... why the limits?

DanPatterson_Retired — Sun, 25 Sep 2016 02:34:22 GMT

Update... as of this date, the geometry class remains the same in ArcMap 10.4.1 and ArcGIS Pro 1.3.1

Maybe some time in the future....

Re: Point class in arcpy... why the limits?

JoshuaBixby — Mon, 26 Sep 2016 14:42:36 GMT

I hear you, frustrating but not surprising. I think the only big improvement/update to the ArcPy Geometry classes that users can expect in the next release is indexing and slicing of geometries.

From conversations I have had, I get the impression there is too much history behind ArcPy Geometry classes to see any major "fixes." Since the Geometry classes have worked a certain way for years (even if that way is less than ideal), and Esri doesn't know all of what their customers do with the classes as is, there is a reluctance to change existing behavior. I can't say I agree completely with the logic when it comes to some of the "fixes," but I have resigned myself, sadly.

I have been told to focus my efforts on the new ArcGIS Python API. Since the API is being written mostly from the ground up, they can (whether they will or not is a different matter) design different behavior than the ArcPy Geometry classes. I can't say I have had much luck, or feel much confidence yet, but I am focusing on trying to convince Esri to aim high with the geometry model in the new API.

Re: Point class in arcpy... why the limits?

DanPatterson_Retired — Mon, 26 Sep 2016 16:50:25 GMT

THe python api unfortunately doesn't address geometry unless the underlying 'classes' have been written. It appears the focus is now on mapping and collecting and not on analysis. I had hopes that when the 3D offerings were enhanced, that would have included the pillars from which they were built

Re: Point class in arcpy... why the limits?

JoshuaBixby — Mon, 26 Sep 2016 17:01:52 GMT

I agree with you on how the API looks now, but I have heard from various folks that analysis (including geometry) will be a future development. I wish I had more or something concrete to offer, but this is one of those times with Esri where I choose to believe because accepting the opposite will be too defeating.