2

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:

enter image description here

  • 2
    Doubles (and floats) are imprecise by their nature. Learn to deal with a very small (some would say insignificant) difference and round it out. If you are expecting to get whole numbers then return int(p.area), if you're expecting small numbers then round to a few decimals. The difference is 0.000244140774% of the total area in this case. – Michael Stimson Mar 16 '15 at 3:00
  • 1
    This may be a coordinate reference issue. What is your default coordinate reference? What is the xy scale of the data frame? There's nothing "forced" about double precision conversion -- that is the correct datatype for coordinate data. – Vince Mar 16 '15 at 3:13
  • @Vince this virtual shape, raster extent expressed in terms of column and row number. no coordinate system – FelixIP Mar 16 '15 at 3:40
  • 1
    You are putting them in as integers but they're stored internally as doubles... The Esri Geometry has to accommodate many different storage systems and has multiple co-classes and interfaces; a compromise has to be made somewhere. In the end no matter how big (wide) you make a floating point it's not going to store all numbers with 100% precision. If you want to calculate the area accurately you'll need to write and implement your own area calculation algorithm or round to an appropriate number of significant digits. It's not that 100 can't be stored as a double, it's the area calculation.. – Michael Stimson Mar 16 '15 at 4:13
  • 3
    @FelixIP Every shape has a coordinate reference. This is intrinsic to ArcGIS. Failure to set a proper coordinate reference (which is much more than a coordinate system) is causing your precision (not accuracy) issue. – Vince Mar 16 '15 at 4:43
0

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
  • Tick, though you answered it already. Although calculating area in GIS using subroutine is way over the top. – FelixIP Apr 21 '15 at 20:15
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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.