Errors in geometry... the importance of Spatial Reference

2828
0
11-17-2015 03:53 AM

Errors in geometry... the importance of Spatial Reference

This thread Errors in arcpy's Polygon class originated a long time ago.

Recently threads have been appearing where the importance of having a spatial reference specified seems to have been forgotten.

The cause is known... the solution is easy, but the lesson isn't sticking.  Read and see for yourself ... as a refresher... or for those venturing into the world of coordinate geometry.

The fix... uncomment line 19  and fix the ending to line 26.  But at least give making a mistake a chance.

geometry_errors_demo.py
start
"""
Script:  geometry_errors_demo.py
Author:  
  Dan.Patterson AT carleton.ca
Modified: 2016-01-30
Purpose:  Demonstrates calculating various properties for polygons
Notes:
  Revised formatting from the original thread:
  
 geonet.esri.com/thread/10256
  
"""
import numpy as np
import arcpy
from textwrap import dedent
triangle = [[0.,0.],[0.,1.],[1.,0.]] 
square = [[0.,0.],[0.,1.],[1.,1.],[1.,0.]] 
rectangle = [[0.,0.],[0.,1.],[2.,1.],[2.,0.]] 
polygons = [triangle, square, rectangle] 
labels = ["Triangle", "Square", "Rectangle"]
#SR = arcpy.SpatialReference(3395)  # uncomment 
# ---- arcpy section ------
def arcpy_demo(polygons, labels):
    """Demo of arcpy errors """
    polys = []
    c = 0
    for p in polygons:
        poly = arcpy.Polygon(arcpy.Array([arcpy.Point(*pair) for pair in p])) #... p]), SR)
        frmt = """
        {} Input coordinates: {}
        Polygon properties...
        Returned coordinates in WKT
        -  {}
        area          {:<20.16f}
        length        {:<20.16f}
        centroid      {:}
        true centroid {:}  
        first point   {:}
        last point    {:}  
        extent -      (XMin, YMin, XMax, YMax, ZMin, ZMax, MMin, MMax)
        - values:     {!s:}
        convex hull:  (x,y), (x,y)...(x,y)
        - sequence:   {:}
        """
        extent_str = str(poly.extent).replace(" ",", ")
        args = [labels[c], polygons[c], poly.WKT, poly.area, poly.length,
                poly.centroid, poly.trueCentroid,
                poly.firstPoint, poly.lastPoint,
                extent_str,poly.hullRectangle]
        print(dedent(frmt).format(*args))
        array.removeAll()
        polys.append(poly)
        c += 1
# ---- part B area determination ----
def ring_area(x, y):
    """x and y as separate lists or arrays from references above
    - assumes first and last points are the same so uses slices
    - exterior rings ordered clockwise
    - holes are ordered counter-clockwise   
    """
    a0 = np.dot(x[1:], y[:-1])
    a1 = np.dot(x[:-1], y[1:])
    area = (a0-a1) * 0.5 # np.abs(a0-a1) * 0.5 for unsigned area, ie no holes
    return area
def poly_area(Xs,Ys):
    """calculate areas for parts"""
    area = 0.0
    for i in range(len(Xs)):
        x = Xs[i]; y = Ys[i]
        area += ring_area(x, y) 
    return area
def area_demo(polygons, labels):
    """Sample area calculations for simple geometry"""
    polygons = [np.array(i) for i in polygons] # convert lists to arrays
    for i in range(len(polygons)):
        x = polygons[i][:,0]; y = polygons[i][:,1]
        area = ring_area(x, y)
        frmt = "\nArea for {}\n{}\nArea... {}"
        print(frmt.format(labels[i],polygons[i], area))
if __name__=="__main__":
    """ """
    arcpy_demo(polygons, labels)
    area_demo(polygons, labels)
    
    

Results
geometry_errors_demo.py outputs
Triangle Input coordinates: [[0.0, 0.0], [0.0, 1.0], [1.0, 0.0]]
Polygon properties...
Returned coordinates in WKT
-  MULTIPOLYGON (((0 0,
                   1.0001220703125 0,
                   0 1.0001220703125,
                   0 0)))
area          0.5001220777630806  
length        3.4146303364895956  
centroid      0.3333740234375 0.3333740234375 NaN NaN
true centroid 0.3333740234375 0.3333740234375 NaN NaN  
first point   0 0 NaN NaN
last point    0 0 NaN NaN  
extent -      (XMin, YMin, XMax, YMax, ZMin, ZMax, MMin, MMax)
- values:     0, 0, 1.0001220703125, 1.0001220703125, NaN, NaN, NaN, NaN
convex hull:  (x,y), (x,y)...(x,y)
- sequence:   1.0001220703125 2.22044604925031E-16
              0.50006103515625 -0.50006103515625
              -0.50006103515625 0.50006103515625
              0 1.0001220703125

Square Input coordinates: [[0.0, 0.0], [0.0, 1.0], [1.0, 1.0], [1.0, 0.0]]
Polygon properties...
Returned coordinates in WKT
-  MULTIPOLYGON (((0 0,
                   1.0001220703125 0, 
                   1.0001220703125 1.0001220703125,
                   0 1.0001220703125,
                   0 0)))
area          1.0002441555261612  
length        4.0004882812500000  
centroid      0.50006103515625 0.50006103515625 NaN NaN
true centroid 0.50006103515625 0.50006103515625 NaN NaN  
first point   0 0 NaN NaN
last point    0 0 NaN NaN  
extent -      (XMin, YMin, XMax, YMax, ZMin, ZMax, MMin, MMax)
- values:     0, 0, 1.0001220703125, 1.0001220703125, NaN, NaN, NaN, NaN
convex hull:  (x,y), (x,y)...(x,y)
- sequence:   6.12398146082414E-17 1.0001220703125 
              1.0001220703125 1.0001220703125 
              1.0001220703125 0 
              0 0

Rectangle Input coordinates: [[0.0, 0.0], [0.0, 1.0], [2.0, 1.0], [2.0, 0.0]]
Polygon properties...
Returned coordinates in WKT
-  MULTIPOLYGON (((0 0,
                   2.0001220703125 0,
                   2.0001220703125 1.0001220703125,
                   0 1.0001220703125,
                   0 0)))
area          2.0003662258386612  
length        6.0004882812500000  
centroid      1.00006103515625 0.50006103515625 NaN NaN
true centroid 1.00006103515625 0.50006103515625 NaN NaN  
first point   0 0 NaN NaN
last point    0 0 NaN NaN  
extent -      (XMin, YMin, XMax, YMax, ZMin, ZMax, MMin, MMax)
- values:     0, 0, 2.0001220703125, 1.0001220703125, NaN, NaN, NaN, NaN
convex hull:  (x,y), (x,y)...(x,y)
- sequence:   6.12398146082414E-17 1.0001220703125 
              2.0001220703125 1.0001220703125 
              2.0001220703125 -2.22044604925031E-16
              0 0

Area for Triangle
[[ 0.  0.]
 [ 0.  1.]
 [ 1.  0.]]
Area... 0.5
Area for Square
[[ 0.  0.]
 [ 0.  1.]
 [ 1.  1.]
 [ 1.  0.]]
Area... 1.0
Area for Rectangle
[[ 0.  0.]
 [ 0.  1.]
 [ 2.  1.]
 [ 2.  0.]]
Area... 2.0

Calculating the area old-school (67-83) proved that sometimes the old methods still have merit.

Version history
Last update:
‎12-12-2021 03:46 AM
Updated by: