Rounded/Erroneous Lat/Long point feature locations using arcpy.da.FeatureClassToNumPyArray

[AllisonBailey]

I created a NumPy array from a CSV file and sued arcpy.da.NumPyArrayToFeatureClass to convert it to a point feature class using the lat/long columns (decimal degrees).   The output location coordinates appear to be rounded, so that many of the features are now co-located.     I checked the input NumPy array values, and they are correct.

Below is a snippet of the first 10 line of the input file (lat/long in 5th/6th columns):

core002,1,6/2/2014,17:41:16,48.563868,-122.922793,-6.018824,-6.018824,0,1,0,0,0
core002,1,6/2/2014,17:41:17,48.563868,-122.922793,-6.248483,-6.248483,0,1,0,0,0
core002,1,6/2/2014,17:41:18,48.563868,-122.922792,-6.346908,-6.346908,0,1,0,0,0
core002,1,6/2/2014,17:41:19,48.563868,-122.922788,-6.478141,-6.478141,0,1,0,0,0
core002,1,6/2/2014,17:41:20,48.563868,-122.922787,-6.839035,-6.839035,0,1,0,0,0
core002,1,6/2/2014,17:41:21,48.563870,-122.922785,-7.167119,-7.167119,0,1,0,0,0
core002,1,6/2/2014,17:41:22,48.563870,-122.922783,-9999.000000,-7.380373,0,1,0,0,0
core002,1,6/2/2014,17:41:23,48.563870,-122.922782,-7.593627,-7.593627,0,1,0,0,0
core002,1,6/2/2014,17:41:24,48.563870,-122.922778,-7.757669,-7.757669,0,1,0,0,0
core002,1,6/2/2014,17:41:25,48.563870,-122.922777,-7.888903,-7.888903,0,1,0,0,0

I compared that to an output NumPy array exported using arcpy.da.FeatureClassToNumPyArray from the point feature class that I created.   All of the XY locations for these 10 records are now identical (and incorrect).   They do change as you go further into the file, but this sample exhibits the co-location issue.

core002,1,6/2/2014,17:41:16,48.563904,-122.922791,-6.018824,-6.018824,0,1,0,0,0
core002,1,6/2/2014,17:41:17,48.563904,-122.922791,-6.248483,-6.248483,0,1,0,0,0
core002,1,6/2/2014,17:41:18,48.563904,-122.922791,-6.346908,-6.346908,0,1,0,0,0
core002,1,6/2/2014,17:41:19,48.563904,-122.922791,-6.478141,-6.478141,0,1,0,0,0
core002,1,6/2/2014,17:41:20,48.563904,-122.922791,-6.839035,-6.839035,0,1,0,0,0
core002,1,6/2/2014,17:41:21,48.563904,-122.922791,-7.167119,-7.167119,0,1,0,0,0
core002,1,6/2/2014,17:41:22,48.563904,-122.922791,-9999.000000,-7.380373,0,1,0,0,0
core002,1,6/2/2014,17:41:23,48.563904,-122.922791,-7.593627,-7.593627,0,1,0,0,0
core002,1,6/2/2014,17:41:24,48.563904,-122.922791,-7.757669,-7.757669,0,1,0,0,0
core002,1,6/2/2014,17:41:25,48.563904,-122.922791,-7.888903,-7.888903,0,1,0,0,0

Anyone have any idea why this is happening.  And, if so, is there a solution?

[JoshuaBixby · Accepted solution]

Using your first 5 coordinates, and assuming WGS 84, I am not seeing any rounding:

>>> import arcpy
>>> import numpy
>>> 
>>> fc = 'in_memory/points'
>>> array = numpy.array([(48.563868,-122.922793),
...                      (48.563868,-122.922793),
...                      (48.563868,-122.922792),
...                      (48.563868,-122.922788),
...                      (48.563868,-122.922787)],
...                      numpy.dtype([('X','<f8'),('Y','<f8')]))
...                      
>>> SR = arcpy.SpatialReference(4326)
>>> arcpy.da.NumPyArrayToFeatureClass(array, fc, ['X','Y'], SR)
>>> desc = arcpy.Describe(fc)
>>> desc.spatialreference.XYResolution
1e-09
>>> desc.spatialreference.XYTolerance
8.983152841195215e-09
>>> arcpy.da.FeatureClassToNumPyArray(fc,"*")
array([(1, [48.56386800000007, -122.92279299999996]),
       (2, [48.56386800000007, -122.92279299999996]),
       (3, [48.56386800000007, -122.92279199999996]),
       (4, [48.56386800000007, -122.92278799999997]),
       (5, [48.56386800000007, -122.92278699999997])], 
      dtype=[('OBJECTID', '<i4'), ('Shape', '<f8', (2,))])
>>> ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

What version of ArcMap or ArcGIS Pro are you using?  What spatial reference are you using?  What is the backend data store for the feature class?

Can you share the code you are using to create the feature class from the NumPy array?

UPDATE:  If you are using Web Mercator (3857), the XY resolution and tolerance will lead to noticeable rounding for your coordinates:

>>> import arcpy
>>> import numpy
>>> 
>>> fc = 'in_memory/points'
>>> array = numpy.array([(48.563868,-122.922793),
...                      (48.563868,-122.922793),
...                      (48.563868,-122.922792),
...                      (48.563868,-122.922788),
...                      (48.563868,-122.922787)],
...                      numpy.dtype([('X','<f8'),('Y','<f8')]))
...                      
>>> SR = arcpy.SpatialReference(3857)
>>> arcpy.da.NumPyArrayToFeatureClass(array, fc, ['X','Y'], SR)
>>> desc = arcpy.Describe(fc)
>>> desc.spatialreference.XYResolution
0.0001
>>> desc.spatialreference.XYTolerance
0.001
>>> arcpy.da.FeatureClassToNumPyArray(fc,"*")
array([(1, [48.563900001347065, -122.92280000075698]),
       (2, [48.563900001347065, -122.92280000075698]),
       (3, [48.563900001347065, -122.92280000075698]),
       (4, [48.563900001347065, -122.92280000075698]),
       (5, [48.563900001347065, -122.92280000075698])], 
      dtype=[('OBJECTID', '<i4'), ('Shape', '<f8', (2,))])
>>> ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[AllisonBailey]

Thanks Joshua.   You nailed it --- it was the spatial reference.    

In my example, I was not assigning any spatial reference, since it is optional.    I guess even without the spatial reference assigned, the XYTolerance and XYResolution are affecting the output feature class.

When I assigned it WGS-84, the coordinates did not get rounded as you pointed out.

[JoshuaBixby]

Although some tools don't require a spatial reference, I would go as far as to say it is a best practice to provide one if the tool allows for or requires it.  There are numerous examples of tools that don't require a spatial reference but behave oddly when one isn't provided. 

[VinceAngelo]

It's best to consider spatial reference use to be as "optional" as breathing.  In reality, every shape has a default spatial reference, which assumes units of meters, and coordinate precision of millimeters.

- V

[AllisonBailey]

Thanks all.   It totally agree that using a spatial reference is a best practice.  The script was a work-in-progress to learn about pandas, numpy and conversion to a feature class, so I wasn't yet concerned with the spatial reference specifics.   Now I have learned, as Joshua demonstrated, that there may be unexpected results if spatial reference is not provided.   And, Vince, I did not know, until now (after ~20 years of experience with ESRI products) that the default ("unknown") spatial reference assumes units of meters and coordinate precision of millimeters.    

[VinceAngelo]

Actually, I should clarify that the units are assumed to be meters (the actual units are undefined) with a tenth-millimeter (1.0e-04) coordinate precision (but what's an order of magnitude between friends?)

The SgShape library became the Esri reference geometry library 20+ years ago, but wasn't immediately adopted by all products. SgShape's CoordRef default origin and scale are {0.0,0.0,1.0},  which would snap degrees to the nearest degree and invalidate anything outside Quadrant I.  The arcpy SpatialReference type manifests the SgShape library's CoordRef as a wrapper, but it sets the nominal origin and scale to {-450359962737.05,-450359962737.05,10000} (which partitions the coordinate space addressable by a long long integer with {0,0} in the center). The arcpy Geometry type retains the double precision array values, so it's only storing them in a feature class that actually truncates the coordinates to four places (as per the SpatialReference)

>>> import arcpy
>>> pt = arcpy.Point()
>>> pt.X = 1.23456789
>>> pt.Y = 9.87654321
>>> print pt
1.23456789 9.87654321 NaN NaN
>>> pg = arcpy.PointGeometry(pt)
>>> sr = pg.spatialReference
>>> print sr.falseOriginAndUnits
-450359962737.05 -450359962737.05 10000
>>> print pg.WKT
POINT (1.2345678899999999 9.8765432099999995)‍‍‍‍‍‍‍‍‍‍‍‍