how to get area for each part of a multipart polygon using python (ArcGIS 9.3)?

Question

Thanks in advance!

Russel

KimOllivier · Answer

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

MarcinGasior · Answer

I'm afraid it's not so easy. The GetPart() method returns Array object not Geometry. And Array doesn't have Area property.

I didn't find easy way to convert back Array to Geometry without using cursors, so I suggest using GP tool Multipart to Singlepart (Data Management) and then calculate area for each single polygon.
If you want to report which multipart polygons were exploder, look at script example in web help for this tool.

Here's the script which uses this tool, writes singlepart feature class to memory (not to disk) and calculates area:import arcgisscripting
gp = arcgisscripting.create(9.3)

infc = r"Z:\Test.gdb\MultipartPolygon"
singlePartFC = r"in_memory\SinglePart"
gp.MultipartToSinglepart_management(infc,singlePartFC)

dsc = gp.Describe(singlePartFC)
shapefieldname = dsc.ShapeFieldName

rows = gp.SearchCursor(singlePartFC)
row = rows.Next()

while row:
 feat = row.getValue(shapefieldname)
 poly_area = feat.Area
 print poly_area

row = rows.Next()

del row, rows

gpnobbir · Answer

You'll get more details in this topic:

http://webhelp.esri.com/arcgisdesktop/9.3/index.cfm?TopicName=Reading_geometries

Here is my modified version of the sample in the above topic:import arcgisscripting
gp = arcgisscripting.create(9.3)

infc = "path to input features"

dsc = gp.Describe(infc)
shapefieldname = dsc.ShapeFieldName

rows = gp.SearchCursor(infc)
row = rows.Next()

while row:

feat = row.getValue(shapefieldname)
    partnum = 0
    partcount = feat.PartCount

# iterate through each part of a multipart feature
    while partnum < partcount:
        poly = feat.GetPart(partnum)
        poly_area = poly.Area
        
        partnum += 1

row = rows.Next()

del row
del rows

The code is not tested thoroughly - you may need to edit if you find any error.