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

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)

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

# specify a named ellipsoid
geod = Geod(ellps="WGS84")

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

area = abs(geod.geometry_area_perimeter(poly)[0])

print('# Geodesic area: {:.3f} m^2'.format(area))

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

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