4

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
3
  • 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.
    – user2856
    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.
    – user2856
    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
    – Big_Mac
    Commented May 15, 2021 at 15:04

2 Answers 2

7

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:

enter image description here

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).

3
  • 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?)
    – Big_Mac
    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.
    – Big_Mac
    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.
    – xunilk
    Commented May 16, 2021 at 2:30
4

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'

#NAD83 / CONUS Albers (https://epsg.io/5070)
crs_conus = 'EPSG:5070'
#NAD83 / Alaska Albers (https://epsg.io/3338)
crs_alaska = 'EPSG:3338'
#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'

gdf = gpd.read_file(shp)

# if we're just looking at CONUS _OR_ Alaska _OR_ Hawaii in one shapefile 
# otherwise need more info about your data to filter and apply different CRSs (i.e. does it have a "STATE" column...)
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_alaska = gdf[(gdf['X'] < -128) & (gdf['Y'] > 49)].to_crs(crs_alaska)
gdf_hawaii = gdf[(gdf['X'] < -128) & (gdf['Y'] < 30)].to_crs(crs_hawaii)

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

new_gdf = gdf_conus.append(gdf_alaska, ignore_index=True).append(gdf_hawaii, ignore_index=True)
new_gdf.drop(['geometry', 'X', 'Y'], axis=1).to_csv('/tmp/test.csv', index=False)

Alternatively if you want to stick with the low level shapely operations, you do not need to explode your polygons to points to reproject:

polygons = [your, list, of, shapely, polygons, for, CONUS]
wgs84 = CRS('EPSG:4326')  # I'm assuming WGS84, but you didn't say (look in the .prj file)
albers_conus = CRS('EPSG:5070')

transformer = Transformer.from_crs(wgs84, albers_conus, always_xy=True).transform
# transformed = [transform(transformer, geom) for geom in polygons]
areas = [transform(transformer, geom).area for geom in polygons]

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.