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
Solved! Go to Solution.
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:
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.
You seem to be missing the getPart, getLength and getArea methods in Polygon class for arcpy
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
Okay how to calculate area for each polygon ?
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.
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:
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.
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.
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...