What a challenge to an old school surveyor's assistant!An array of coordinate pairs is what we always had after reducing a traverse, when we had to calculate the parcel area.This was hard to do with only an adding machine and 10 figure log book. Hence the Double Difference Area method: Consider each segment is a trapezoid above the X axis.The area is half times the average of the Y of each side times the difference of the X values.Sum these areas and that is the net area.Wait until the end to halve the Y values once.As you go around the polygon, coming in reverse has negative values under the polygon.I'm not quite sure why Esri gets a different result (1 in a million), must be the round off of the precision.# area of multiparts
import arcgisscripting
gp = arcgisscripting.create(9.3)
infc = "C:/path.gdb/multipolyfc"
dsc = gp.Describe(infc)
shapefieldname = dsc.ShapeFieldName
rows = gp.SearchCursor(infc) #,"multipart = 2")
row = rows.Next()
while row:
feat = row.getValue(shapefieldname)
partnum = 0
partcount = feat.PartCount
partArea = 0
cumArea = 0
print partcount ,feat.area
# iterate through each part of a multipart feature
while partnum < partcount:
arrayPoly = feat.GetPart(partnum)
# iterate through each vertex of the array
# use old pre-computer Double Difference Area survey trick
pt = arrayPoly.next()
x0 = pt.x
y0 = pt.y
pt = arrayPoly.next()
partArea = 0
while pt:
partArea += (pt.x - x0) * (pt.y + y0)
x0 = pt.x
y0 = pt.y
pt = arrayPoly.next()
partnum += 1
partArea = partArea / 2.0
print partnum,partArea
cumArea += partArea
print "count %d,cum %f, feat %f diff (%12.8f) " % (partcount, cumArea,feat.area ,(feat.area - cumArea)/100)
# break
row = rows.Next()
del row
del rows