10

I'm serializing my arcpy geometries as geojson so that I can 'hydrate' them back as geometries later and I'm having 2 problems in the cycle.:

PROBLEM 1: Precision

    R0 = arcpy.SearchCursor(self.shpTest, "FID=0").next().getValue("Shape")          
    geojson = R0.__geo_interface__                        
    R1 = arcpy.AsShape(geojson)
    self.assertTrue(R0.equals(R1)) <<< THIS FAILS

If I check the string representation, the coordinates have slightly changed:

    geojson2 = R1.__geo_interface__
    print geojson
    print geojson2  

    {'type': 'Polygon', 'coordinates': [[(442343.5516410945, 4814166.6184399202), (442772.17749834526, 4811610.7383281607), (441565.67508534156, 4811499.6131059099), (440772.50052100699, 4814184.7808806188), (442343.5516410945, 4814166.6184399202)]]}
    {'type': 'Polygon', 'coordinates': [[(442343.55169677734, 4814166.6185302734), (442772.17749023438, 4811610.73828125), (441565.67510986328, 4811499.6130981445), (440772.50048828125, 4814184.7808837891), (442343.55169677734, 4814166.6185302734)]]}

PROBLEM 2: Holes If the polygon has holes, geo_interface generates an error:

    R0_WithHoles = arcpy.SearchCursor(self.shpTest, "FID=0").next().getValue("Shape")          
    geojson = R0.__geo_interface__  <<< generates this ERROR:

    File "C:\Program Files\ArcGIS\Desktop10.0\arcpy\arcpy\arcobjects\geometries.py", line 68, in __geo_interface__
        return {'type': 'Polygon', 'coordinates': [[(pt.X, pt.Y) for pt in part] for part in self]}
    AttributeError: 'NoneType' object has no attribute 'X'

Any ideas on how to solve these problems?

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
Víctor Velarde
  • 383
  • 1
  • 13
  • Yep, just ran across number 2 myself. And doesn't appear to be much love for this topic. – Random Geo guy Dec 08 '11 at 00:43
  • This is still broken in arcpy in ArcGIS 10.1 -- It would be nice if ESRI could comment on the subject. – James Mills Sep 17 '12 at 23:42
  • I came across the first and second problems. With me, the coordinated do not seem to change (when you print them) but geom1.equals(geom2) fails me only a few times. I am not sure why that happens as well. The 2nd problem was fixed using @valveLondon 's suggestion. If you found out how to fix the .equals please do share. – Michalis Avraam Nov 07 '12 at 17:46
  • @MichalisAvraam We had the same issue also and got onto ESRI for a solution - turns out it's a known bug (when you create a geom without a projection it truncates the precision) - take a look at this question also. – om_henners Nov 12 '12 at 21:12
  • @om_henners I assumed that. But the arcpy.AsShape() function does not let you specify a spatial reference. I have set all environment variables hoping it would do something (output coords, etc...). The solution then is to manually decode the GeoJSON because ESRI doesn't care about accuracy? – Michalis Avraam Nov 12 '12 at 23:40
  • @MichalisAvraam Afraid so. Not sure about 10.1, but ESRI have said they're not going to fix it for 10. If it helps I've got a script that I've just put up on gist.github that you might be able to strip for relavent bits. – om_henners Nov 13 '12 at 00:01
  • @om_henners Thank you. I wound up rewriting the AsShape function in geometries.py, but ended up not using it and switching to ogr to handle geometries. But this is incredibly helpful. Thank you so much for sharing. Hopefully others can see this too. – Michalis Avraam Nov 13 '12 at 19:51

1 Answers1

5

OK - well I thought I had solved it.

replace line ~80 of this file C:\Python26\ArcGIS10.0\Lib\arcpy\arcobjects\geometries.py from this:

return {'type': 'Polygon', 'coordinates': [[(pt.X, pt.Y) for pt in part] for part in self]}

to this (or something that is more concise and elegant and does the same thing):

  obj = {"type": "Polygon"}
    coordinates = []
    for part in self:
        _part = []
        for pt in part:
            if pt is not None:
                print pt
                _part.append([pt.X,pt.Y])
            else:
                print "none"
                coordinates.append(_part)
                _part=[]
        coordinates.append(_part)
    obj["coordinates"]=coordinates
    return obj

Basically they forgot to consider donuts in the shape which are marked by null point values. This spits out good geoJson (separate parts) but the arcpy.AsShape method trashes GeoJSON.

this code:

import arcpy
gj = {
  'type': 'Polygon', 'coordinates': [
   [[-122.803764, 45.509158], [-122.796246, 45.500050], [-122.808193, 45.500109],
      [-122.803764, 45.509158]],
   [[-122.804206, 45.504509], [-122.802882, 45.502522], [-122.801866, 45.504479], 
      [-122.804206, 45.504509]]
   ]
 }

 p = arcpy.AsShape(gj)
 print p.__geo_interface__

outputs this:

    {'type': 'Polygon', 'coordinates': [[[-122.8037109375, 45.50927734375],  
    [-122.79620361328125, 45.5001220703125], [-122.80810546875, 45.5001220703125],
    [-122.8037109375, 45.50927734375]]]}

I give up. ;)

Update The holes problem has been solved at 10.1 with this chunk of python:

return {'type': 'Polygon', 'coordinates': [[((pt.X, pt.Y) if pt else None)
                                                    for pt in part]
                                                        for part in self]}
Random Geo guy
  • 2,010
  • 3
  • 19
  • 31
  • shouldn't that return a dictionary instead of a string representing a dictionary? :) – blah238 Mar 13 '12 at 01:14
  • Yes, you're right, it should. I changed it to spit out a valid GeoJSON dictionary obj. but after checking against the AsShape method I realized the futility of my efforts. – Random Geo guy Mar 14 '12 at 22:58
  • I wonder if that has anything to do with the problem described in this thread: http://forums.arcgis.com/threads/9763-Errors-in-arcpy-s-Polygon-class -- should have been fixed in 10 SP2 and definitely 10.1. – blah238 Mar 14 '12 at 23:36
  • Yah - that looks like a pretty solid thread. I just put SP4 on and .. no love. – Random Geo guy Mar 14 '12 at 23:51
  • @valveLondon- Did Esri update their master code for geometries.py in 10.1 so that it now works "out of the box"? Or does the user still have to change the chunk of code on their own? – RyanKDalton Jun 25 '12 at 21:33
  • 2
    ESRI updated C:\Program Files\ArcGIS\Server\arcpy\arcpy\arcobjects\geometries.py at 10.1 but if you're at 10.0 you can fix it yourself. – Random Geo guy Jun 25 '12 at 21:35
  • 3
    Yes, I fixed it in 10.1, the update above is the new source in the .py file. I thought it made it into a service pack for 10 but I guess not. – Jason Scheirer Jun 25 '12 at 22:03
  • Still hope for SP5, chin up Jason! – blah238 Jun 25 '12 at 23:51
  • @valveLondon Unfortunately that breaks the Polygon class in 10.0 (the _fromGeoJson gives an error when called from arcpy.AsShape() because it cannot handle None types). Do you mind copy/pasting that part too? The return statement. Now it is return cls(Array([map(lambda p: Point(*p), part) for part in coordinates])) – Michalis Avraam Nov 06 '12 at 22:10
  • @valveLondon My solution follows: return cls(Array([map(lambda p: Point(*p) if p is not None else Point(), part) for part in coordinates])). I hope it is as efficient as the potential ESRI solution. – Michalis Avraam Nov 06 '12 at 23:27
  • @JasonScheirer- Do you know if the fix made it into 10SP5? And if the change I have in my comment would solve the AsShape part of the geometries.py file? – Michalis Avraam Nov 08 '12 at 22:25