# Calculate area and perimeter of coordinate points in Python

I have a set of points that make up a polygon. I need to find the area and perimeter for the shape in Python. I have used the `shapely` library to calculate the area and perimeter but it is nonsensical because the points are latitude and longitude coordinates (in US States).

I need to convert the points to an equal-area projection system first. I was thinking of using Sinusoidal projection. I was able to do this and calculate area (code below), but I'm not sure how to interpret the result (also is it right?). What unit is it in? (m²?)

Next I need to calculate perimeter, but would the Sinusoidal projection distort it and throw it off? How would I deal with that?

I am a computer scientist and my work on a project has led me down this rabbit hole. How can I calculate the area/perimeter of this shape in the correct way?

``````def sinusoidal(lat, long):
earth_radius = 6371009 # in meters
lat_dist = pi * earth_radius / 180.0

y = lat * lat_dist
x = long * lat_dist * cos(radians(lat))

return x, y
``````
• You've barely even looked into the entrance to the rabbit hole. Actually, it's more like a bottomless pit :) How accurate do you need your area calcs to be? That transformation assumes the earth is a sphere and it's not, there are better projections that can be used. If you're concerned about accuracy, I suggest you look at `pyproj` if you want to stick with low level shapely or geopandas which is a higher level library. Commented May 14, 2021 at 23:17
• If you post some example data (do you just have points or have you extracted them from something else that could be directly read by a geospatial library) we might be able to suggest another way. Commented May 14, 2021 at 23:18
• @user2856 I have a shape file of about 500 shapes in the following formats: DBF, PRJ, SHP, SHX Commented May 15, 2021 at 15:04

In this link you have an example for determining Geodesic area for a polygon located in antarctic region by using pyproj python3 module. Adapting such script for following polygon (green) in an arbitrary area of USA:

it looks as follows:

``````from pyproj import Geod

geod = Geod('+a=6378137 +f=0.0033528106647475126')

lons = [-102.20253523916332, -101.59096157206567, -100.65438018473898, -101.90199046561818]
lats = [37.21550522238942, 37.70825886273666, 36.93243398218993, 36.394249155143996 ]

poly_area, poly_perimeter = geod.polygon_area_perimeter(lons, lats)

print("area: {} , perimeter: {}".format(poly_area, poly_perimeter))
``````

After running it in python console, result was as follow:

``````area: 10473941945.624039 , perimeter: 418057.54140898824
``````

Obtained values are comparable (differences probably due to different used Geod in each case) to those showed in attributes table of green shapefile (area: 10476912480.580376, perimeter: 418057.541405565).

• Thanks so much, this is super helpful! Is there a different Geod you would recommend for Hawaii/Alaska or should the increased error not really matter much? Also what units is the perimeter and area being returned in? (meters/meters^2?) Commented May 15, 2021 at 21:34
• Deeper dive into the documentation mentions it is in fact in meters/meters^2. Still curious about Hawaii in particular. Commented May 15, 2021 at 21:41
• You're welcome. After importing Geod, you can do help(Geod). There is an extensive documentation with many example usage (including shapely use). Additional important information can also be found with a help(pyproj); after importing pyproj. Commented May 16, 2021 at 2:30

The following simple geopandas script (a higher level library that makes use of shapely, fiona/gdal & pandas) shows how to reproject to an Albers Equal Area CRS and calculate the area (in the units of the CRS which in this case is metres):

``````import geopandas as gpd

shp = '/tmp/test.shp'

crs_conus = 'EPSG:5070'
#NAD83 / Hawaii Albers Equal Area Conic (https://epsg.io/102007)  *not an EPSG crs
crs_hawaii = '+proj=aea +lat_1=8 +lat_2=18 +lat_0=13 +lon_0=-157 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs'

# if we're just looking at CONUS _OR_ Alaska _OR_ Hawaii in one shapefile
gdf = gdf.to_crs(crs_conus) # or crs_alaska etc

gdf['area_m'] = gdf['geometry'].area
gdf.drop(['geometry'], axis=1).to_csv('/tmp/test.csv', index=False)
``````

Or if CONUS, AK and HI in one shapefile with no attributes to specify State

``````# a very simple spatial filter

gdf['X'] = gdf['geometry'].centroid.x
gdf['Y'] = gdf['geometry'].centroid.y

gdf_conus = gdf[gdf['X'] > -128].to_crs(crs_conus)
gdf_hawaii = gdf[(gdf['X'] < -128) & (gdf['Y'] < 30)].to_crs(crs_hawaii)

gdf_conus['area_m'] = gdf_conus['geometry'].area
gdf_hawaii['area_m'] = gdf_hawaii['geometry'].area

``````polygons = [your, list, of, shapely, polygons, for, CONUS]