Calculate Polygon area in planar units (e.g. square-meters) in Shapely

Question

I am using Python 3.4 and shapely 1.3.2 to create a Polygon object out of a list of long/lat coordinate pairs which I transform into a well-known-text string in order to parse them. Such a Polygon might look like:

POLYGON ((-116.904 43.371, -116.823 43.389, -116.895 43.407, -116.908 43.375, -116.904 43.371))

Since shapely does not handle any projections and implements all geometry objects in the carthesian space, calling the area method on that polygon like:

poly.area

gives me the area of that polygon in the unit of square-degrees. To get the area in a planar unit like square-meters, I guess that I would have to transform the coordinates of the polygon using a different projection (which one?).

I read several times that the pyproj library should provide the way to do this. Using pyproj, is there a way to transform a whole shapely Polygon object into another projection and then calculate the area?

I do some other stuff with my Polygons (not what you think now) and only in certain cases, I need to calculate the area.

So far, I only found this example: http://all-geo.org/volcan01010/2012/11/change-coordinates-with-pyproj/

which would mean splitting each Polygon object into its outer and, if present, inner rings, grab the coordinates, transform each pair of coordinates into another projection and rebuild the Polygon object, then calculate its area (what unit is it then anyway?). This looks like a solution, but is not very practical.

Any better ideas?

Stomachic

Answer 1

Calculate a geodesic area, which is very accurate and only requires an ellipsoid (not a projection). This can be done with pyproj 2.3.0 or later.

from pyproj import Geod
from shapely import wkt

poly = wkt.loads("""\
    POLYGON ((-116.904 43.371, -116.823 43.389, -116.895 43.407,
              -116.908 43.375, -116.904 43.371))
""")

print(f"Cartesian area: {poly.area:.6f}")
# Cartesian area: 0.001467

# Specify a named ellipsoid
geod = Geod(ellps="WGS84")
geod_area = abs(geod.geometry_area_perimeter(poly)[0])

print(f"Geodesic area: {geod_area:.3f} m^2")
# Geodesic area: 13205034.647 m^2

abs() is used to return only positive areas. A negative area may be returned depending on the winding direction of the polygon.

Chasitychasm

Nice. Aside, from the pyproj documentation: For polygons, holes should use the opposite traversal to the exterior (if the exterior is CCW, the holes/interiors should be CW). – Adaptable

Very helpful. As a shapely newb I was left wondering what was up with .area returning 0. – Simferopol

Answer 2

Ok, I finally made it with the Basemap toolkit of the matplotlib library. I will explain how it works, maybe this will be helpful for somebody sometime.

1. Download and install the matplotlib library on your system. http://matplotlib.org/downloads.html

For Windows binaries, I recommend using this page: http://www.lfd.uci.edu/~gohlke/pythonlibs/#matplotlib Beware of the hint that says:

Requires numpy, dateutil, pytz, pyparsing, six, and optionally pillow, pycairo, tornado, wxpython, pyside, pyqt, ghostscript, miktex, ffmpeg, mencoder, avconv, or imagemagick.

Hence, if not already installed on your system, you have to download and install numpy, dateutil, pytz, pyparsing and six as well for matplotlib to work properly (for Windows users: all of them can be found on the page, for Linux users, the python package manager "pip" should do the trick).

2. Download and install the "basemap" toolkit for matplotlib. Either from http://matplotlib.org/basemap/users/installing.html or for Windows binaries also from here: http://www.lfd.uci.edu/~gohlke/pythonlibs/#basemap

3. Do the projection in your Python code:

#Import necessary libraries
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np

#Coordinates that are to be transformed
long = -112.367
lat = 41.013

#Create a basemap for your projection. Which one you use is up to you.
#Some examples can be found at http://matplotlib.org/basemap/users/mapsetup.html
m = Basemap(projection='sinu',lon_0=0,resolution='c')

#Project the coordinates:
projected_coordinates = m(long,lat)

Output:

projected_coordinates (10587117.191355567, 14567974.064658936)

Simple as that. Now, when using the projected coordinates to build a polygon with shapely and then calculating the area via shapely's area method, you'll get the area in the unit of square-meters (according to the projection you used). To get square-kilometers, divide by 10^6.

Edit: I struggled hard not to transform only single coordinates, but whole geometry objects like Polygons when those were already given as shapely objects and not via their raw coordinates. This meant writing a lot of code to

• get the coordinates of the outer ring of the polygon
• determine if the polygon has holes and if so, process each hole seperately
• transform each pair of coordinates of the outer ring and any holes
• put the whole thing back together and create a polygon object with the projected coordinates
• and that's only for polygons... add even more loops to this for multipolygons and geometrycollections

Then I stumbled upon this part of the shapely documentation which makes things a lot easier: http://toblerity.org/shapely/manual.html#shapely.ops.transform

When the projection map is set, for example as done above:

m = Basemap(width=1,height=1, resolution='l',projection='laea',lat_ts=50,lat_0=50,lon_0=-107.)

Then, one can transform any shapely geometry object using this projection via:

from shapely.ops import transform
projected_geometry = transform(m,geometry_object)

Stomachic

See also example for projection transformation here. Note that the accuracy of the area calculation depends on the projection. E.g., use a UTM zone if the data all fits nicely in one. – Chasitychasm

Answer 3

Convert to radians and assuming the Earth is perfect sphere of 6370Km radius:

p = np.array([[-116.904,43.371],
              [-116.823, 43.389],
              [-116.895,43.407],
              [-116.908,43.375],
              [-116.904,43.371]])

poly = Polygon(np.radians(p))

poly.area
# 4.468737548271707e-07

poly.area*6370**2
# 18.132751662246623

Johannisberger

This one is very wrong, it assumes that the polygon is very small & on the equator. Othwerwise, the math is way more complicated than this! (@DST, I suggest you delete or improve this post, as it is massively misleading.) – Peng