In reply to Robert.

Yes, with badly scaled data, it is best to do the computation in a transformed coordinate system in which the coordinates have unit scale and then transform it back. If you want a more principled approach to replace your ad-hoc calculations, try using this rewritten module for the polygon centroid:

/* _PolyCentroid: Return the centroid of a simple polygon.
P is an N x 2 matrix of (x,y) values for consectutive vertices. N > 2. */
start _PolyCentroid(P);
lagIdx = 2:nrow(P) || 1; /* vertices 2, 3, ..., N, 1 */
xi = P[ ,1]; yi = P[ ,2]; /* x_i and y_i for i=1..N */
/* transform to standardized coordinates, do the computation, then restore */
meanX = mean(xi); meanY = mean(yi); /* mean of data */
stdX = std(xi); stdY = std(yi); /* std dev of data */
xi = (xi-meanX)/stdX; yi = (yi-meanY)/stdY; /* standardize */

xip1 = xi[lagIdx]; yip1 = yi[lagIdx]; /* x_{i+1} and y_{i+1}, x_{N+1}=x_1 */
d = xi#yip1 - xip1#yi; /* vector of difference of products */
A = 0.5 * sum(d); /* signed area of polygon */
cx = sum((xi + xip1) # d) / (6*A); /* X coord of centroid */
cy = sum((yi + yip1) # d) / (6*A); /* Y coord of centroid */
/* convert back to original coordinates */
cx = stdX*cx+meanX; cy = stdY*cy+meanY;
return( cx || cy );
finish;