The issue is most likely one of floating point precision. In your case you've already extracted the geometries using arcpy, and you've matched them with your RUID.
Happily since you've got arcpy installed you've got numpy, which makes comparing sets of numeric arrays easy. In this case I'd suggest the numpy.allclose function, which is available in numpy 1.3.0 (installed with ArcGIS 10).
From the samples you gave above
geom1geoJSON = {'type': 'Polygon', 'coordinates': [[(-122.8423481559999, 47.060497293000083), (-122.84239755599992, 47.059262423000064), (-122.84416913599989, 47.059309693000046), (-122.84416913599989, 47.060497293000083), (-122.8423481559999, 47.060497293000083)]]}
geom2geoJSON = {'type': 'Polygon', 'coordinates': [[(-122.8423481559999, 47.060497293000083), (-122.84239755599992, 47.059262423000064), (-122.84416913599989, 47.059309693000046), (-122.84416913599989, 47.060497293000083), (-122.8423481559999, 47.060497293000083)]]}
import numpy as np
close = np.allclose(np.array(geom1geoJSON["coordinates"]), np.array(geom2geoJSON["coordinates"]), atol=1e-7)
#Returns True
The atol
keyword specifies the tolerance value.
Note that you shouldn't use arcpy.AsShape
at all. Ever. As I noted in this question (/shameless plug) there's a known bug in ArcGIS that truncates geometries when they are created without a coordinate system (even after setting the env.XYTolerance
environment variable). In arcpy.AsShape
there's no way to avoid this. Fortunately geometry.__geo_interface__
does extract the correct geometries from an existing geometries (though it doesn't handle complex polygons without the fix from @JasonScheirer ).
I did this using a python script and numpy to help me with the matrix algebra for rotation and translation of every point.
Here is the key function. You have to unpack the vertices from the feature, transform and reassemble the feature.
import arcpy
import math
import numpy
def trans(px,py,tx,ty,angle):
'''
Rotate a point px,py around origin 0,0
Sense clockwise or anticlockwise
(default clockwise)
Transform to plot coords tx,ty
Does not handle z or m values
'''
# angle in degrees +/- 360
th = math.radians(angle) # theta (radians)
if sense == 'anticlock':
matrix = [[math.cos(th),-math.sin(th),tx],
[ math.sin(th), math.cos(th),ty],
[0,0,1]]
elif sense == 'clock':
matrix = [[ math.cos(th),math.sin(th),tx],
[-math.sin(th),math.cos(th),ty],
[0,0,1]]
else:
raise Exception, "sense error in DrawPlotR"
a = numpy.array(matrix)
b = numpy.array([px,py,1])
c = numpy.dot(a,b)
return (c[0],c[1])
Best Answer
This code should do it using the SHAPE@XY token that came with arcpy.da.UpdateCursor in ArcGIS 10.1.
The coding pattern used here came from ArcPy Café.