None isn't...nor is 0 or 1 ... more explorations into geometry

DanPatterson_Retired · ‎04-17-2015

I will leave this as a muse to which I will return. I thought I would post it now so I can return to more important issues regarding standardization and geometry expression. Writing this chapter is going to be a pain...

Anyway...here goes....

Dan

When null/nothing/nadda/zilch is not that.

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

Oh well, it has to be better elsewhere.

>>> null_point.__geo_interface__
Traceback (most recent call last):
  File "<interactive input>", line 1, in <module>
AttributeError: 'Point' object has no attribute '__geo_interface__'
>>> 
>>> null_geometry = arcpy.PointGeometry(null_point)
>>> null_geometry.__geo_interface__
{'type': 'Point', 'coordinates': (0.0, 0.0)}
>>>
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Nope...'errors' at worst, misrepresentation at best. I will try the reverse later.

So I began to experiment with other formats before readdressing this issue to see
whether WKT, JSON offered any alternatives. So...I foolishly thought, lets explore
some real shape before returning to this issue. The foolish journey begins

>>> import arcpy
>>> SR = arcpy.SpatialReference(4326)  # alway use a coordinate system ... GCS_WGS_1984
>>> pnts = [[0.0,0.0],[0.0,1.0],[1.0,0.0],[0.0,0.0]]  # make a polygon ... a triangle
>>> arc_poly = arcpy.Polygon(arcpy.Array([arcpy.Point(*xy) for xy in pnts]), SR) # now a Polygon
>>> # ... now for some magic ...
>>> arc_poly.__geo_interface__ # check the Polygon class for more info (sort of)
{'type': 'MultiPolygon', 'coordinates': [[[(5.684341886080802e-14, 5.684341886080802e-14), (5.684341886080802e-14, 1.0000000000000568), (1.0000000000000568, 5.684341886080802e-14), (5.684341886080802e-14, 5.684341886080802e-14)]]]}
>>> # ?????? 
>>> Esri_dict = {"rings": [pnts], "spatialReference" : {"wkid" : 4326}}
>>> Esri_json = arcpy.AsShape(Esri_dict,True)
>>> Esri_json.JSON
u'{"rings":[[[0,0],[0,1],[1,0],[0,0]]],"spatialReference":{"wkid":4326,"latestWkid":4326}}'
>>> WKT_poly = arcpy.FromWKT(pnts_str,SR)
>>> WKT_poly.WKT
u'MULTIPOLYGON (((0 0, 1 0, 1 1, 0 1, 0 0)))'
>>> # where did the SR go???
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Sigh... there is more to explore in the example script below. Standardization is just

'''
formats_to_polygon.py
Author:  Dan.Patterson@carleton.ca

see:
https://geonet.esri.com/thread/122755

https://geonet.esri.com/thread/144517

https://geonet.esri.com/thread/63090

http://en.wikipedia.org/wiki/Well-known_text

http://gis.stackexchange.com/questions/48949/

           epsg-3857-or-4326-for-googlemaps-openstreetmap-and-leafletothers:
Notes:
- http://resources.arcgis.com/en/help/main/10.2/index.html#//03q300000066000000

    "AsShape does not support GeoJSON coordinate reference system objects; the
    spatial reference of a geometry object created from GeoJSON will be Unknown."
- arcpy.Polygon.__geo_interface__  see...
https://gist.github.com/sgillies/2217756
 and https://github.com/cleder/pygeoif

http://w3facility.org/question/arcpy-geometry-__geo_interface
__
         -and-asshape-function-loss-of-precision-and-holes/
Coordinate systems:
- 4326: http://spatialreference.org/ref/epsg/wgs-84/

        Horizontal component of 3D system. Used by the GPS satellite navigation
        system and for NATO military geodetic surveying.
        Google Earth
- 4326 = GEOGCS["GCS_WGS_1984",
                DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],
                PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]
- 3857:  http://spatialreference.org/ref/sr-org/7483/
  or SR-ORG:7483
        Google Maps
- 3857 = PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["GCS_WGS_1984",
                DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],
                PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],
                PROJECTION["Mercator"],PARAMETER["central_meridian",0],
                PARAMETER["scale_factor",1],PARAMETER["false_easting",0],
                PARAMETER["false_northing",0],UNIT["Meter",1]]
- 3395: http://spatialreference.org/ref/epsg/3395/

         Very small scale mapping. Last Revised: June 2, 2006,
         World - between 80°S and 84°N
- 3395 = PROJCS["WGS 84 / World Mercator",
                GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",
                SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],
                UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator"],
                PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],
                PARAMETER["false_easting",0],PARAMETER["false_northing",0],
                UNIT["Meter",1]] 
'''
import arcpy
import json
def main(code):
  pnts = [[0.0,0.0],[0.0,1.0],[1.0,0.0],[0.0,0.0]]
  pnts_str = "POLYGON((0.0 0.0, 1.0 0.0, 0.0 1.0, 0.0 0.0))"
  Esri_dict = {"rings": [pnts], "spatialReference" : {"wkid" : code}}
  Geo_dict = {"type": "Polygon","coordinates": [pnts],
                  "spatialReference":{"wkid":code,"latestWkid":code}}
  SR = arcpy.SpatialReference(code)
  arc_poly = arcpy.Polygon(arcpy.Array([arcpy.Point(*xy) for xy in pnts]), SR) 
  WKT_poly = arcpy.FromWKT(pnts_str,SR)
  Esri_json = arcpy.AsShape(Esri_dict,True)
  Geo_json = arcpy.AsShape(Geo_dict,False)
  print('\nTriangle Input formats for .... {}, {}'.format(SR.name,SR.factoryCode))
  print(' arcpy.Polygon:\n {}'.format(repr(pnts)))       
  print(' WKT:\n  input:  {}\n  output: {}  '.format(pnts_str,WKT_poly.WKT))
  print(' Esri JSON:\n  input:  {}\n  output: {}  '.format(Esri_dict,Esri_json.JSON))
  print(' GeoJSON:\n  input:  {}\n  output:   {}  '.format(Geo_dict,Geo_json.JSON))
  print(' Geo_json.__geo_interface__ :\n  output:  {} '.format(Geo_json.__geo_interface__))
#-------------------------------------------------------------------------
if __name__=='__main__':
  for code in [3395,3857,4326]:  # modify the list to suit your purposes
    main(code)
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Enjoy... as I said, there is way more to this, but this post is too long already.

JoshuaBixby · ‎04-18-2015

I don't mean to get hung up on semantics, but I usually can't help myself. Are you interested in a null point or empty point? Without getting into the debate of whether Python has a null or not, isn't the Python equivalent of a null anything None?

Even when a company like Esri, Microsoft, Oracle, or [insert open source product here] says they are implementing the OGC Simple Feature standards (Simple Feature Access - Part 1: Common Architecture, Simple Feature Access - Part 2: SQL Option, etc...), I have come to expect idiosyncrasies when it comes to how those standards are implemented or translate to a company's geometry data type.

When it comes to Esri, specifically ArcPy, it took me a while to figure out that they do in fact implement an empty point or POINT EMPTY. My initial confusion came from not appreciating what I was reading in Esri's documentation: Point -- Help | ArcGIS for Desktop.

Discussion
A Point is not a geometry class, but is commonly used to construct geometry. In the example below, a Point is used to create a PointGeometry object.

It is a bit odd that Esri puts the Point documentation in with the geometry classes, but the documentation itself does clearly state that a Point is not a geometry class. Typing the objects in Python does affirm what the documentation states:

>>> point = arcpy.Point()
>>> type(point)
<class 'arcpy.arcobjects.arcobjects.Point'>
>>> pointGeometry = arcpy.PointGeometry(point)
>>> type(pointGeometry)
<class 'arcpy.arcobjects.geometries.PointGeometry'>

It is at this point that I have to start assuming what Esri was thinking because I have yet to find any documentation that explains the following.

Since Point is not a geometry class, it cannot implement an empty point, or at least it doesn't. Instead, an empty point is implemented in PointGeometry. Where things are really awkward is that PointGeometry as a point geometry constructor requires either Point or Array objects. If Point can't construct an empty point, how does one go about making PointGeometry construct an empty point? Maybe pass an empty Array object or an Array object with a None object?

>>> point
<Point (0.0, 0.0, #, #)>
>>> pointGeometry.WKT
u'POINT (0 0)'
>>> pointGeometry = arcpy.PointGeometry(arcpy.Array())
Runtime error 
Traceback (most recent call last):
    ....
RuntimeError: Object: CreateObject cannot create geometry from inputs
>>> pointGeometry = arcpy.PointGeometry(arcpy.Array(None))
Runtime error 
Traceback (most recent call last):
    ....
RuntimeError: Object: CreateObject cannot create geometry from inputs

Nope, that doesn't work. Fortunately, there are other geometry constructors. My preference is to use arcpy.FromWKT (FromWKT -- Help | ArcGIS for Desktop).

>>> pointGeometry = arcpy.FromWKT('POINT EMPTY')
>>> pointGeometry.WKT
u'POINT EMPTY'
>>> pointGeometry.JSON
u'{"x":"NaN","y":"NaN","spatialReference":{"wkid":null}}'
>>> repr(pointGeometry.firstPoint)
'None'

If the Point associated with an empty point is really of NoneType instead of arcpy.arcobjects.arcobjects.Point, can we pass a None object to PointGeometry to create an empty point?

>>> arcpy.PointGeometry(None)
Runtime error 
Traceback (most recent call last):
    ....
RuntimeError: Object: CreateObject cannot create geometry from inputs

Again, no. The question now, or at least when I first reached this point, is whether ArcGIS and ArcPy actually implement an empty point or whether the empty point from using arcpy.FromWKT just looks like an empty point. I can't say I have done exhausting testing, but a few cursory checks indicate it isn't just appearances.

>>> point_empty = arcpy.FromWKT('POINT EMPTY')
>>> point_empty.WKT
u'POINT EMPTY'
>>> point = arcpy.FromWKT('POINT (1 1)')
>>> point.WKT
u'POINT (1 1)'
>>> #Check that union of point with empty point is point
>>> point.union(point_empty).WKT
u'MULTIPOINT ((1 1))'
>>> #Check that intersect of point with empty point is empty point
>>> point.intersect(point_empty, 1).WKT
u'POINT EMPTY'
>>> #Check that point contains empty point
>>> point.contains(point_empty)
True
>>> #Check that point is disjoint from empty point
>>> point.disjoint(point_empty)
True

In the end, it seems ArcGIS and ArcPy implement an empty point, but the ArcPy Geometry Classes fail as constructors to support it.

DanPatterson_Retired · ‎04-18-2015

Sorry...can't get syntax highlighting on for some reason today...

Joshua...

Basically we are on the same page, whether you realize it or not. Point...PointGeometry...Whatever.

This is fundamentally wrong

>>> import arcpy

>>> arcpy.Point()

>>> print arcpy.PointGeometry(arcpy.Point()).WKT

POINT (0 0)

>>>

since 0,0 for coordinates, regardless of their implementation, are valid observations.

rolling out my own is more fun especially if you can trap arcpy.Point() and substitute it with..

>>> my_point = arcpy.Point(np.NAN,np.NAN)
>>> exactly = my_point.X
>>> exactly
nan
>>> type(exactly)
<type 'float'>
>>> exactly == 0.0
False
>>> exactly < 0.0
False
>>> exactly > 0.0
False
>>> exactly == np.Inf
False
>>> exactly < np.Inf
False
>>> exactly > np.Inf
False
>>> exactly == None
False
>>> exactly != None
True
>>>

Fixed format

JoshuaBixby · ‎04-18-2015

I agree. I had some thoughts in my head on this subject, and your blog post was the nudge I needed to get them written down.

Whether a bug or design failure, I think something is remiss with the ArcPy Geometry Classes when it comes to constructing geometries. If Esri wants to claim that Point is not a geometry class so it can't support an empty point, fine, but then arcpy.Point() and arcpy.Point(None) should throw exceptions instead of returning POINT (0 0). Additionally, PointGeometry should support calling with no arguments or a None argument to generate an empty point instead of throwing exceptions the way it does now.

Just for fun, it is interesting to see how Shapely handles this situation:

>>> from shapely.geometry import Point
>>> point_empty = Point()
>>> type(point_empty)
<class 'shapely.geometry.point.Point'>
>>> point_empty.wkt
'GEOMETRYCOLLECTION EMPTY'

With Shapely, the Point object is actually a geometry class, unlike with ArcPy. Also, calling Point with no arguments returns an empty point as an empty geometry collection.

JoshuaBixby · ‎04-18-2015

Dan Patterson, regarding syntax highlighting in blog comments, it is a hassle because the Advanced Editor isn't available. I typically compose my blog comments in a draft discussion post and then copy and paste the raw HTML into the blog comment. Annoying but gets the job done.

DanPatterson_Retired · ‎04-18-2015

Ahhhh I will give that a shot...as soon as it says I am the owner of the blog...grief