# Calculating area in km² for Polygon in WKT using Python

I have polygons in WKT that looks like this

``````POLYGON ((24.8085317 46.8512821, 24.7986952 46.8574619, 24.8088238 46.8664741, 24.8155239 46.8576335, 24.8085317 46.8512821))
``````

I would like to calculate the area of this polygon in km².

How do I do this in Python?

It wasn't readily apparent to me how to use @sgillies answer, so here is a more verbose version:

``````import pyproj
import shapely
import shapely.ops as ops
from shapely.geometry.polygon import Polygon
from functools import partial

geom = Polygon([(0, 0), (0, 10), (10, 10), (10, 0), (0, 0)])
geom_area = ops.transform(
partial(
pyproj.transform,
pyproj.Proj(init='EPSG:4326'),
pyproj.Proj(
proj='aea',
lat_1=geom.bounds,
lat_2=geom.bounds
)
),
geom)

# Print the area in m^2
print(geom_area.area)
``````

It looks like your coordinates are longitude and latitude.

Use Shapely's `shapely.ops.transform` function to transform the polygon to projected equal area coordinates and then take the area.

``````python
import pyproj
from functools import partial

geom_aea = transform(
partial(
pyproj.transform,
pyproj.Proj(init='EPSG:4326'),
pyproj.Proj(
proj='aea',
lat1=geom.bounds,
lat2=geom.bounds)),
geom)

print(geom_aea.area)
# Output in m^2: 1083461.9234313113
``````

The other answers seem to be correct, except that at some point, the lat1 and lat2 parameters in the pyproj code were renamed with underscores: lat_1 and lat_2 (see https://stackoverflow.com/a/55259718/1538758).

``````import pyproj
import shapely
import shapely.ops as ops
from shapely.geometry.polygon import Polygon
from functools import partial

geom = Polygon([(0, 0), (0, 10), (10, 10), (10, 0), (0, 0)])
geom_area = ops.transform(
partial(
pyproj.transform,
pyproj.Proj(init='EPSG:4326'),
pyproj.Proj(
proj='aea',
lat_1=geom.bounds,
lat_2=geom.bounds)),
geom)

# Print the area in m^2
print geom_area.area
``````

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

POLYGON ((24.8085317 46.8512821, 24.7986952 46.8574619,
24.8088238 46.8664741, 24.8155239 46.8576335, 24.8085317 46.8512821))''')

area = abs(geod.geometry_area_perimeter(poly))

print('# Geodesic area: {:12.3f} m²'.format(area))
print('#                {:12.3f} km²'.format(area/1e6))

# Geodesic area:  1083466.869 m²
#                       1.083 km²
``````

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

I've stumbled across "area" which seems simpler to use:

https://pypi.org/project/area/

For example:

``````from area import area

obj = {'type':'Polygon','coordinates':[[[24.8085317,46.8512821], [24.7986952,46.8574619], [24.8088238,46.8664741], [24.8155239,46.8576335], [24.8085317,46.8512821]]]}

area_m2 = area(obj)

area_km2 = area_m2 / 1e+6
print ('area m2:' + str(area_m2))
print ('area km2:' + str(area_km2))
``````

... returns:

area m2:1082979.880942425

area km2:1.082979880942425

I've found that the projected option is a bit too slow for my use so I did some digging and found a javascript lib from Mapbox called cheap-ruler. The "cheap" in the name is intended to imply that this algorithm isn't intended to return the right answer, but simply a good enough one.

Below is my copy and paste port of the area function to work with a shapely geom. My port doesn't remove area for holes, but it could be updated to do so.

``````import math

def area(geom):
## Setup from:
RE = 6378.137
FE = 1 / 298.257223563
E2 = FE * (2 - FE)

lat = geom.centroid.y
m = RAD * RE * 1000

w2 = 1 / (1 - E2 * (1 - coslat * coslat))
w = math.sqrt(w2)
kx = m * w * coslat
ky = m * w * w2 * (1 - E2)

## Area calc from:
ring = geom.exterior.coords

sumVal = 0
j = 0
l = len(ring)
k = l - 1
while j < l:
sumVal += wrap(ring[j] - ring[k]) * (ring[j] + ring[k])
k = j
j += 1

return (abs(sumVal) / 2) * kx * ky;

def wrap(deg):
while deg < -180:
deg += 360
while deg > 180:
deg -= 360
return deg
``````

As you probably use this in a loop for many geometries you can get faster transformations by using modern Proj as intended and by not recreating static objects like the initial CRS:

``````from pyproj import Transformer, transform, CRS
from shapely.geometry import Polygon
import shapely.ops as ops

geom = Polygon([(0, 0), (0, 10), (10, 10), (10, 0), (0, 0)])

crs_4326 = CRS.from_epsg(4326)

transformer = Transformer.from_crs(
crs_4326,
CRS(proj='aea',
lat_1=geom.bounds,
lat_2=geom.bounds
)
)

geom_area = ops.transform(transformer.transform, geom)
``````

is about twice as fast for me if you measure the setup of the transformer and the transformation.

@bugmenot123 answer is not correct you need to change code to otherwise you get wrong area you will have switched long,lat ....

``````        transformer = Transformer.from_crs(
CRS_WGS_84,
CRS(proj='aea',
lat_1=shape.bounds,
lat_2=shape.bounds
),
always_xy=True
)
``````