nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

Problem in calculation Area and Perimeter using Python

AhsanAbbas

I need a script that determines the perimeter (in meters) and area (in square meters) of each of the individual islands of the Hawaii.shp feature class. This shapefile is multipart feature class.

My code illustrate below, but it doesn't give me correct answer ?

fc = r"D:\Python Exercise\Exercise08-2014-12-05\Exercise08\Hawaii.shp"

cursor = arcpy.da.SearchCursor(fc, ["OID@","SHAPE@","SHAPE@AREA","SHAPE@LENGTH"])

units = arcpy.Describe(fc).spatialReference.linearUnitName

for row in cursor:

print "Total Area: ", row[2]

print("Feature {0}: ".format(row[0]))

partnum = 0

perimeter = 0

for part in row[1]:

print("Part {0}: ".format(partnum))

for point in part:

perimeter += row[3]

partnum += 1

#print "Area: ", row[2]

print "Perimeter: ", perimeter, units

Hawaii.zip

Find more posts tagged with

Accepted answers

JoshuaBixby

Given the number and nature of your questions today, not to mention the references to "exercise" in path names, it seems you are relying pretty heavily on GeoNet users to figure your class exercises out for you. I don't mind helping a new GIS user get started with arcpy and Python, but I don't get the sense the Help documentation has been thoroughly read before asking for others to provide code.

That said, and if I understand your question correctly, I believe the following will get close to what you are after:

fc = #path for feature class
SR = arcpy.Describe(fc).spatialReference
units = SR.linearUnitName

cursor = arcpy.da.SearchCursor(fc, ["OID@","SHAPE@"])
for OID, shape in cursor:
    print "Feature: {}\tTotal Area: {} {}^2".format(
        OID, shape.area, units
    )
    partnum = 0
    for part in shape:
        pn = arcpy.Polygon(part, SR)
        print "\tPart: {}\t Area: {} {}^2\tPerimeter: {} {}".format(
            partnum, pn.area, units, pn.length, units
        )
        partnum += 1

A few comments to move beyond simply giving you a proverbial fish:

If you are going to retrieve the shape as a geometry object, there isn't much point in using the other shape tokens. Since retrieving a shape as a geometry object carries the most overhead, the other tokens were created to speed up the search cursor when one didn't need the full shape but only its area, length, etc....
The code above, since it relies on the area properties instead of the getArea method of the polygon object, assumes you are interested in planar values and using a feature class that is already projected. If you are interested in other measurement types or are using unprojected data, then you will want to look at the getArea method.
When retrieving the parts of a polygon, what is returned aren't polygons but arrays of points that need to be reconstructed back into polygons if you are interested in measuring areas and perimeters.

I strongly encourage you to read the Polygon (arcpy) help, I think it will help clarify your misunderstanding of how to work with arcpy geometry objects.

All comments

AhsanAbbas

Actually i didn't mean to write code for me, i provided my code above, i just want to make it correct, where i am doing wrong nothing else...

JoshuaBixby

Another tool to consider is Multipart To Singlepart (Data Management).‌ If you took this approach, then you could use the SHAPE@AREA and SHAPE@LENGTH tokens and not get involved directly with manipulating arcpy geometry objects.

JoshuaBixby

That said, and if I understand your question correctly, I believe the following will get close to what you are after:

fc = #path for feature class
SR = arcpy.Describe(fc).spatialReference
units = SR.linearUnitName

cursor = arcpy.da.SearchCursor(fc, ["OID@","SHAPE@"])
for OID, shape in cursor:
    print "Feature: {}\tTotal Area: {} {}^2".format(
        OID, shape.area, units
    )
    partnum = 0
    for part in shape:
        pn = arcpy.Polygon(part, SR)
        print "\tPart: {}\t Area: {} {}^2\tPerimeter: {} {}".format(
            partnum, pn.area, units, pn.length, units
        )
        partnum += 1

A few comments to move beyond simply giving you a proverbial fish:

If you are going to retrieve the shape as a geometry object, there isn't much point in using the other shape tokens. Since retrieving a shape as a geometry object carries the most overhead, the other tokens were created to speed up the search cursor when one didn't need the full shape but only its area, length, etc....
The code above, since it relies on the area properties instead of the getArea method of the polygon object, assumes you are interested in planar values and using a feature class that is already projected. If you are interested in other measurement types or are using unprojected data, then you will want to look at the getArea method.
When retrieving the parts of a polygon, what is returned aren't polygons but arrays of points that need to be reconstructed back into polygons if you are interested in measuring areas and perimeters.

I strongly encourage you to read the Polygon (arcpy) help, I think it will help clarify your misunderstanding of how to work with arcpy geometry objects.

RichardFairhurst

Well I misread your earlier post. !Shape@Area! and !Shape@Length! are the area and perimeter values of the feature in its native units. So is this just a simple unit conversion problem? If so, once you know the native units of the feature you can get the conversion factors from here and apply them to the row[2] and row[3] values to get square meters and meters. Or you could provide a Spatial Reference that used meters as the cursor's fourth parameter and then the !Shape@Area! and !Shape@Length! would return values in square meters and meters.

If you are trying to independently calculate area and confirm the values against those of the area and length returned by the shape method, than the euclidean distance method I gave will be one method for doing that for the perimeter length, provided the distances are not covering a large portion of the globe's surface. For area calculations I don't know an independent formula that will work for all possible polygon shapes (which in part accounts for why we buy Esri's product to take advantage of people who have figured that out). Using the !Shape@Area! values, or the getArea values from a geometry would be the approaches I would take.

AhsanAbbas

Okay how to calculate area for each polygon ?

RichardFairhurst

You appear to be cumulatively adding the total shape length of the entire polygon to the perimeter variable for each point along the perimeter that you encounter. The length between points would require a euclidean distance calculation between all point pairs that make up the part being added as you traverse the part.

The Shape.Length field (row[3]) should already be the perimeter of the polygon if this shapefile was converted from a geodatabase feature class or calculated using the geometry calculator. Did you mean to update the Shape.Length field after determining the perimeter? If so you need to use an UpdateCursor, not a SearchCursor.

I believe something more like this untested code should work if you wanted to get the euclidean distances from a shapefile that is using a Spatial Reference with linear units (however, Dan's option of using getLength from the geometry may be all you need to do):

for row in cursor:
  # Create the geometry object
  #
  feat = row[1]

  partnum = 0
  perimeter = 0
  # Step through each part of the feature
  #
  for part in feat:
    # Print the part number
    #
    print "Part %i:" % partnum

    lastpnt = None

    # Step through each vertex in the feature
    #
    for pnt in feat.getPart(partnum):
      if pnt:
        # Print x,y coordinates of current point
        #
        print pnt.X, pnt.Y
        if lastpnt == None:
          lastpnt = pnt
        else:
          perimeter += math.sqrt((pnt.X - lastpnt.X) ** 2 + (pnt.Y - lastpnt.Y) ** 2)
      else:
        # If pnt is None, this represents an interior ring
        #
        print "Interior Ring:"
        lastpnt = None
    partnum += 1

DanPatterson_Retired

You seem to be missing the getPart, getLength and getArea methods in Polygon class for arcpy

Join the Conversation

Signed in members can post, follow updates, and more. New here? Register a free account.

Quick Links

Community Help
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.