Accuracy loss when creating polygons with Polygon()

973
3
Jump to solution
10-08-2014 05:41 AM
ChrisHills
Occasional Contributor

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?

0 Kudos
1 Solution

Accepted Solutions
DanPatterson_Retired
MVP Emeritus

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

View solution in original post

3 Replies
DanPatterson_Retired
MVP Emeritus

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

ChrisHills
Occasional Contributor

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

0 Kudos
ChrisHills
Occasional Contributor

Thank you Dan for your helpful answer.

0 Kudos