The Shape.ReturnCenter request
returns the centroid of a shape's extent, not of the shape
itself. To see this, simply create a rectangle, cut it diagonally, and compute
the "Center" of either of the triangular halves. You will find that the
center lies along the hypotenuse, not in the interior where it belongs.
For many applications you need a bona fide centroid, or center of
mass. Mark Cederholm kindly contributed the "return true centroid of
polygon" script to the ESRI web page many months ago (your solution #2),
observing at the time that ReturnCenter does not generally produce a centroid.
However, that script has several problems, including failure to
recognizeexceptional cases (zero-area and null polygons) and a limitation
(properly noted) to one-part polygons. I thought to stand on Mark's shoulders,
as it were, and offer a small modification of his script to fix these
problems; but somehow I screwed up his formula in so doing. Therefore I offer a
solution using a different formula--at least as simple as Mark's--that
I have verified to work. This solution can be incorporated in any of
the scripts you referenced simply by replacing expressions like
[shape].ReturnCenter with expressions like av.Run("Polygon.ReturnCentroid",
[shape]).
--Bill Huber
Quantitative Decisions
Merion, PA
===========================================================
'
' Polygon.ReturnCentroid
'
' Despite what ArcView's online help says, Polygon.ReturnCenter does NOT
' return the true centroid. (It appears to return the centroid of the
' polygon's "extent," or bounding box.)
'
' This routine implements the line integrals of 0.5*x^2*dy and -0.5*y^2*dx,
' which by Green's theorem compute the (area) integrals of x*dx*dy and y*dx*dy,
' respectively. These, when divided by the area (computed as the line integral
' of x*dy), produce the coordinates of the centroid, by definition. Such line
' integrals are the sum of integrals over the edges. Parametrizing an edge
' from (x0,y0) to (x1, y1) as t |--(x0, y0) + t*(x1-x0, y1-y0), each edge
' integral becomes the (ordinary) integral from 0 to 1 of the expression
' created by the following substitutions:
' x --x0 + t*(x1-x0)
' y --y0 + t*(y1-y0)
' dx --(x1-x0) * dt
' dy --(y1-y0) * dt
'
' After substitution, each integral becomes a polynomial in t, homogeneous in
' x0,x1,y0,and y1. According to elementary calculus (integral of dt is t, of t*dt is
' t^2/2, etc.), integration effectively replaces dt by 1, t*dt by 1/2, t^2*dt by 1/3, etc.
'
' For example, the area expression is computed through two substitutions:
' (1) x * dy --(x0 + t*(x1-x0)) * (y1-y0) * dt
' (2) --(x0 + 1/2 * (x1-x0) * (y1-y0) = 1/2 * (x0+x1) * (y1-y0).
'
' Since line integrals respect orientation, one simply sums the appropriate
' integrals over each polygon component, then normalizes by signed area at
' the end to get the centroid. This is possible because ArcView consistently
' represents polygons with appropriately oriented boundaries: see the help for
' Polygon.Clean.
'
' William Huber, 1/26/99 (Tested in AV 3.1).
' whuber@quantdec.com
'
' Arguments: thePolygon (NOT {thepolygon} !)
' Returns: ptCentroid
'
' Comments: A point is always returned. It will be a null point whenever
' an invalid, null, or zero-area polygon is offered as the argument.
'-------------------------------------------------------------------------------'
if (SELF.Is(Polygon).Not or SELF.IsNull) then return point.MakeNull
end
lstlstPts = SELF.Clean.AsList
ptC = 0@0 ' Will be six times the area times the centroid
A = 0 ' Will be twice the (signed) area
ptO = SELF.ReturnCenter ' For numerical stability, use this as the origin.
for each lstPts in lstlstPts ' Process each part
if (lstPts.Count < 1) then continue end
ptX0 = lstPts.Get(lstPts.Count-1) - ptO
for each ptX1 in lstPts ' Integrate over each edge ptX0--ptX1
ptX1 = ptX1 - ptO
x0 = ptX0.GetX
y0 = ptX0.GetY
x1 = ptX1.GetX
y1 = ptX1.GetY
dx = x1-x0
dy = y1-y0
sx = x1+x0
sy = y1+y0
A = sx * dy + A
ptC = ((sx^2 - (x0*x1)) * dy) @ ((y0*y1 - (sy^2)) * dx) + ptC
ptX0 = ptX1
end
end
if (A = 0) then
return point.MakeNull
else
A = 3*A
return ptC / (A@A) + ptO
end
' end of script