Cómo manejar la precisión de las coordenadas en ArcGIS

Question

Estoy creando un polígono rectangular muy básico usando un código Python muy básico

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)

y este es un mensaje

10000.0244140774

Parece que la conversión forzada a dobles es una razón, como se describe en Las exportaciones de ArcGIS 10 pierden precisión con los números grandes

¿Quién puede decir cómo obtener un resultado más preciso?

P.D. El polígono creado en ArcView 3 utilizando las mismas coordenadas tiene un área de 10000.0000

Esto es una actualización de mi pregunta original . Código ligeramente modificado:

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)

produjo la respuesta correcta (!). Área del polígono = 10000.

Esto es un alivio, no tengo que escribir mi propio procedimiento en GIS para calcular el área de la forma. Lo que no me gusta son las coordenadas de los puntos que forman este polígono:

Answer 1

Las coordenadas se convierten en dobles cuando el arcpy Polygon y sin un sistema de referencia de coordenadas adecuado, ArcGIS no sabe qué precisión debe tener el resultado. Si se utiliza un sistema de coordenadas adecuado, es probable que el error disminuya; sin embargo, nunca desaparecerá debido a la limitada precisión del formato doble.

Además (sin saber todo sobre tu caso de uso), usar filas y columnas de un raster de entrada como sistema de referencia es cuanto menos sospechoso. Yo sugeriría o bien utilizar una herramienta alternativa para lo que sea que estés tratando de lograr con la imagen no geocodificada (como la Python Imaging Library, la caja de herramientas de procesamiento de imágenes de MATLAB o cualquier SIG basado en raster) o geocodificar la imagen primero y luego realizar la acción como el recorte.

Answer 2

Esta observación se remonta a 2010. Si no se especifica una referencia espacial, se obtendrán resultados de menor precisión. Véase https://geonet.esri.com/thread/10256 para empezar. Los ejemplos son más extensos. Hay aún más en los años siguientes. Especificar una entrada de punto flotante no tiene ningún impacto, así que descarta esa idea.

>>> 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

¿Qué tal el Python puro? (ejemplo muy simple, obviamente)

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))) 

resultado...

>>> Unit square area: 1.0