Skip navigation
All People > Dan_Patterson > Py... blog > 2015 > April
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.