topic Accuracy loss when creating polygons with Polygon() in Python Questions

Accuracy loss when creating polygons with Polygon()

ChrisHills — Wed, 08 Oct 2014 12:41:31 GMT

When creating polygons with Polygon() the results are not accurate as expected. Take the following script:-

import arcpy


startx = 285000
starty = 138000
boxsize = 1536
boxcount = 2

features = []

for colx in range(1,boxcount+1):
  for coly in range(1,boxcount+1):
    top    = (starty+(boxsize*coly))
    bottom = (starty+(boxsize*(coly-1)))
    left   = (startx+(boxsize*(colx-1)))
    right  = (startx+(boxsize*colx))
    feature_info = [[[left, bottom],
                     [right, bottom],
                     [right, top],
                     [left, top]]]
    for feature in feature_info:
      features.append(arcpy.Polygon(arcpy.Array([arcpy.Point(*coords) for coords in feature])))
      print left, bottom, right, top

for feature in features:
  print feature.extent

This produces the following output:-

285000 138000 286536 139536

285000 139536 286536 141072

286536 138000 288072 139536

286536 139536 288072 141072

285000.00012207 138000.00012207 286536.00012207 139536.00012207 NaN NaN NaN NaN

285000.00012207 139536.00012207 286536.00012207 141072.00012207 NaN NaN NaN NaN

286536.00012207 138000.00012207 288072.00012207 139536.00012207 NaN NaN NaN NaN

286536.00012207 139536.00012207 288072.00012207 141072.00012207 NaN NaN NaN NaN

As you can see, instead of the resulting boundary being equal to the input points, it has changed them by 0.00012207. This is unexpected and unwanted. How do I create my polygons without this error?

Re: Accuracy loss when creating polygons with Polygon()

DanPatterson_Retired — Wed, 08 Oct 2014 13:36:29 GMT

That is an old issue... have a read Create a spatial reference (aka projection) for the data. any UTM zone will work

Re: Accuracy loss when creating polygons with Polygon()

ChrisHills — Thu, 09 Oct 2014 06:41:03 GMT

Thank you Dan for your helpful answer.

Re: Accuracy loss when creating polygons with Polygon()

ChrisHills — Thu, 09 Oct 2014 09:24:42 GMT

I updated my script to include a spatial reference and it now gives the expected result:-

import arcpy  
  
  
startx = 285000  
starty = 138000  
boxsize = 1536  
boxcount = 2  
  
features = []  

# British National Grid
srid = arcpy.SpatialReference(27700)

for colx in range(1,boxcount+1):  
  for coly in range(1,boxcount+1):  
    top    = (starty+(boxsize*coly))  
    bottom = (starty+(boxsize*(coly-1)))  
    left   = (startx+(boxsize*(colx-1)))  
    right  = (startx+(boxsize*colx))  
    feature_info = [[[left, bottom],  
                     [right, bottom],  
                     [right, top],  
                     [left, top]]]  
    for feature in feature_info:  
      features.append(arcpy.Polygon(arcpy.Array([arcpy.Point(*coords) for coords in feature]),srid))  
      print left, bottom, right, top  
  
for feature in features:  
  print feature.extent  
285000 138000 286536 139536
285000 139536 286536 141072
286536 138000 288072 139536
286536 139536 288072 141072
285000 138000 286536 139536 NaN NaN NaN NaN
285000 139536 286536 141072 NaN NaN NaN NaN
286536 138000 288072 139536 NaN NaN NaN NaN
286536 139536 288072 141072 NaN NaN NaN NaN