How to handle coordinates accuracy in ArcGIS

Question

I am creating very basic rectangular polygon using very basic a Python code

nCols,nRows=100,100
corners =[[0,0],[0,nRows],[nCols,nRows],[nCols,0],[0,0]]
p=arcpy.Polygon(arcpy.Array([arcpy.Point(*coords) for coords in corners]))
arcpy.AddMessage(p.area)

and this is a message

10000.0244140774

It seems that forced conversion to doubles is a reason, as described in ArcGIS 10 exports losing precision with large numbers

Who can please tell how to get more accurate result?

P.S. Polygon created in ArcView 3 using same coordinates has an area of 10000.0000

This is update on my original question. Slightly modified code:

nCols,nRows= 101, 101
corners =[[1,1],[1,nRows],[nCols,nRows],[nCols,1],[1,1]]
p=[arcpy.PointGeometry(arcpy.Point(*coords)) for coords in corners]
arcpy.CopyFeatures_management(p, "d:/rubbish/points.shp")
p=arcpy.Polygon(arcpy.Array([arcpy.Point(*coords) for coords in corners]))
arcpy.AddMessage(p.area)

produced correct answer(!). Area of polygon = 10000.

This is a relief, I don't have to write my own procedure in GIS to calculate area of the shape. What I don't like are coordinates of the points that makes this polygon:

Answer 1

The coordinates get converted to doubles when the arcpy Polygon object is created and without a proper coordinate reference system, ArcGIS does not know how precise the result should be. If you use an appropriate coordinate system, the error is likely to diminish; however, it will never disappear due to limited double format precision.

Besides (not knowing everything about your use case), using rows and columns of an input raster as the reference system is suspicious at least. I'd suggest either using an alternative tool for whatever you're trying to accomplish with the non-geocoded image (such as the Python Imaging Library, MATLAB Image Processing toolbox or any raster-based GIS) or geocoding the image first and then performing the action such as clipping.

Answer 2

This observation goes back to 2010. If a spatial reference isn't specified, then you will get lower precision results. See https://geonet.esri.com/thread/10256 for starters. The examples are more extensive. There are even more in subsequent years. Specifying a floating point input has no impact, so rule that out as a thought.

>>> import arcpy
>>> corners =[[0.0,0.0],[0.0,1.0],[1.0,1.0],[1.0,0.0],[0.0,0.0]]
>>> p=arcpy.Polygon(arcpy.Array([arcpy.Point(*coords) for coords in corners]))
>>> p.area
1.0002441555261612
>>> # hmmmm close... try again
>>> SR = arcpy.SpatialReference(3395) # WGS_1984_ World_Mercator_3395.prj doesn't really matter
>>> p=arcpy.Polygon(arcpy.Array([arcpy.Point(*coords) for coords in corners]),SR)
>>> p.area
1.0

How about pure Python? (very simple example obviously)

def polygon_area(pnts):  
  '''determine the area given a list of points'''  
  area = 0.0;  n = len(pnts)  
  j = n - 1;  i = 0  
  for point in pnts:  
    p0 = pnts[i];  p1 = pnts[j]  
    area += p0[0] * p1[1]  
    area -= p0[1] * p1[0]  
    j = i; i += 1  
  area /= 2.0  
  return area 

if __name__=='__main__':
  square = [[0,0],[0,1],[1,1],[1,0]]  
  print('Unit square area: {}'.format(polygon_area(square))) 

result...

>>> Unit square area: 1.0