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