ArcGIS Python API's WKT doesn't match EWKT from PostGIS

[DavidAskov1]

Hi all - I am trying to retrieve features from an ArcGIS Server 10.31 map service using the ArcGIS Python API and insert them into PostgreSQL/PostGIS. For the geometry, both use WKT, but they don't have the same format. Here is a snippet of what I run to get the data:

import arcgis
from arcgis.gis import GIS
from arcgis.features import FeatureLayer
from arcgis.geometry import Geometry

perimetersLAYER = FeatureLayer('https://<URL>/arcgis/rest/services/<service>/MapServer/1')
perimeters_feature_set = perimetersLAYER.query(where="objectid = 8605")

for feature in perimeters_feature_set.features:
    print('\n\n')
    for field in feature.fields:
        print(field + ': ' + str(feature.get_value(field)))
    
    geom = Geometry(feature.geometry)
    print('geom.WKT: ' + str(geom.WKT))
    ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

The ArcGIS Python API Geometry response looks like this:

MULTIPOLYGON (((-13645793.9305409509688616,  4645922.6947647286579013), 
(-13645734.9949775058776140,  4645907.8913778141140938), 
(-13645738.5948890540748835,  4645918.5636747097596526), 
(-13645706.8243041206151247,  4645943.8052548961713910), [SNIP]

I'd like to pass this directly to a SQL insert with the PostGIS MULTIPOLYGON function, which looks like this:

MULTIPOLYGON(((0 0 0,4 0 0,4 4 0,0 4 0,0 0 0),(1 1 0,2 1 0,2 2 0,1 2 0,1 1 0)),((-1 -1 0,-1 -2 0,-2 -2 0,-2 -1 0,-1 -1 0)))

Note that the geometry formats are slightly different. In PostGIS, the x/y/z coords are space delimited and each vertex is comma delimited. The Esri geometry's x/y coords are comma delimited, and each vertex is wrapped in parentheses and also comma delimited.

I can parse all the strings and munge the formats from one to another, but is there a way to automagically convert from one to another? It looks like Esri's WKT does not match the input PostGIS needs for its WKT functions. When I try the Geometry's WKB and JSON properties, it returns None, so I'm not sure where to go with that.

Of course the geometry format is the minutia of my problem. The bigger picture issue is to find the shortest path from an ArcGIS Enterprise map service to PostGIS. If there's a better way, I'm all ears, though I'd like to keep it in python and run it on Linux.

thanks!

[JoshuaBixby]

Since EWKT and EWKB are more a PostGIS creation than a standard like WKT and WKB, I would use the latter if you are actually moving data between systems.  I don't have anything against EWKT or EWKB, I just don't see them as very robust from a systems interoperability perspective.

This isn't an overall Esri issue but one for the ArcGIS API for Python.  If you look at ArcPy, WKT is handled properly:

>>> pg = arcpy.FromWKT('MULTIPOLYGON Z(((0 0 0,4 0 0,4 4 0,0 4 0,0 0 0),(1 1 0,2 1 0,2 2 0,1 2 0,1 1 0)),((-1 -1 0,-1 -2 0,-2 -2 0,-2 -1 0,-1 -1 0)))')
>>> pg.WKT
u'MULTIPOLYGON Z (((0 0 0, 4 0 0, 4 4 0, 0 4 0, 0 0 0), (1 1 0, 1 2 0, 2 2 0, 2 1 0, 1 1 0)), ((-1 -1 0, -2 -1 0, -2 -2 0, -1 -2 0, -1 -1 0)))'
>>> ‍‍‍‍‍‍‍‍

The WKT coming from the Python API is clearly butchered, in a couple of ways.  I haven't investigated enough to know whether this is a unique situation or a broader issue in the code.

Sharing with https://community.esri.com/groups/arcgis-python-api?sr=search&searchId=a9b70294-0900-4a50-b48c-ae5a9...‌

[DavidAskov1]

Joshua - Thanks. I'm not too familiar with the delta between WKT & EWKT. For simple features, they seem to be about the same. Setting that issue aside, however, the ArcGIS Python API's WKT is quite different, so thanks for confirming I'm not crazy or doing something wrong in my code.

I appreciate the ArcPy solution, but we're hoping to avoid all that overhead for now on our little Linux server. Without ArcPy, is it safe to say that writing some code to reformat ArcGIS WKT into PostGIS WKT is our fastest approach? It's not a big job, but then I suppose we'll have to eventually do point, line, and polygon, plus the "multi" versions of each. Any other suggestions are welcome.

[JoshuaBixby]

What version of the API are you using?  I created a simple polygon in one of my hosted feature layers, and I am not seeing the odd WKT output.

>>> import arcgis
>>> print(arcgis.__version__)
1.2.4
>>>
>>> from arcgis.gis import GIS
>>> from arcgis.features import FeatureLayer
>>> from arcgis.geometry import Geometry
>>>
>>> perimetersLAYER = FeatureLayer('https://fqdn/arcgis/rest/services/Hosted/SNF_LandStatus/FeatureServer/0')
>>> perimeters_feature_set = perimetersLAYER.query(where='objectid=4')
>>> feature = next(f for f in perimeters_feature_set)
>>> geom = Geometry(feature.geometry)
>>> print(geom.WKT)
MULTIPOLYGON (((463774.0959999999 5387170.5483999997, 474346.17119999975 5213428.8544999994, 751430.00449999981 5214531.3746000007, 733288.65990000032 5391894.7993000001, 513147.3208999997 5387069.2414999995, 463774.0959999999 5387170.5483999997)))
>>>‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[DavidAskov1]

I am also running 1.2.4. I tried it on a Jupyter notebook and a regular install of Python on a Macbook - same result.

The ArcGIS Server is 10.31, could that be the isssue?

[JoshuaBixby]

Unless your organization has a public-facing service, or you can find a public-facing service that demonstrates the issue, I don't have any more specific suggestions.  I believe my organization still has some 10.3.1 services running, so I will test on one of them when I can find the time.

UPDATE:  I found a public-facing 10.3.1 server my organization has with some feature services.  The following feature service outputs normal WKT using ArcGIS API for Python.

>>> from arcgis.gis import GIS
>>> from arcgis.features import FeatureLayer
>>> from arcgis.geometry import Geometry
>>>
>>> url = "https://apps.fs.usda.gov/arcx/rest/services/EDW_FEATURE/EDW_Wilderness_01/FeatureServer/0"
>>> perimetersLAYER = FeatureLayer(url)
>>> perimeters_feature_set = perimetersLAYER.query(where='wid=708')
>>> feature = next(f for f in perimeters_feature_set)
>>> geom = Geometry(feature.geometry)
>>> print(geom.WKT)
>>> geom.WKT
'MULTIPOLYGON (((-80.309682070000008 37.883744509999985, -80.309682090000024 37.883581180000022, ..., -80.309682070000008 37.883744509999985)))'‍‍‍‍‍‍‍‍‍‍‍‍

[DavidAskov1]

How strange. I am just getting the same non-WKT I posted in my original question in both a Jupyter notebook and on a MacBook running python 3.6. Thanks for your help, but for now, we just dumped it to JSON, read that into a Python dictionary, and transmogrified it into WKT that PostGIS will accept.

I may circle back on this in Windows using the ArcGIS-supplied python just out of curiosity, but if I can't count on it working reliably in all environments (this will probably go to an Ubuntu Server or AWS Lambda for production), then I can't use it.