28

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?

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
amaatouq
  • 441
  • 1
  • 6
  • 10

8 Answers8

30

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[1],
            lat_2=geom.bounds[3]
        )
    ),
    geom)


Print the area in m^2


print(geom_area.area)
Community
  • 1
jczaplew
  • 1,208
  • 12
  • 19
20

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[1],
        lat2=geom.bounds[3])),
geom)


print(geom_aea.area)


Output in m^2: 1083461.9234313113
PolyGeo
  • 65,136
  • 29
  • 109
  • 338
sgillies
  • 9,056
  • 1
  • 33
  • 41
10

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[1],
            lat_2=geom.bounds[3])),
    geom)


Print the area in m^2


print geom_area.area
PolyGeo
  • 65,136
  • 29
  • 109
  • 338
MartinThurn
  • 203
  • 2
  • 5
8

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 ((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)[0])


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.

Taras
  • 32,823
  • 4
  • 66
  • 137
Mike T
  • 42,095
  • 10
  • 126
  • 187
6

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

Mike Honey
  • 613
  • 4
  • 10
2

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:
    ## https://github.com/mapbox/cheap-ruler/blob/48ad4768a52dc176b01494d090cce19f02c7afdd/index.js#L71-L82
    RE = 6378.137
    RAD = math.pi / 180
    FE = 1 / 298.257223563
    E2 = FE * (2 - FE)


lat = geom.centroid.y
m = RAD * RE * 1000
coslat = math.cos(lat * RAD);

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

## Area calc from:
## https://github.com/mapbox/cheap-ruler/blob/48ad4768a52dc176b01494d090cce19f02c7afdd/index.js#L185-L197
ring = geom.exterior.coords

sumVal = 0
j = 0
l = len(ring)
k = l - 1
while j < l:
    sumVal += wrap(ring[j][0] - ring[k][0]) * (ring[j][1] + ring[k][1])
    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
Matthew Whittemore
  • 121
  • 2
2

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[1],
        lat_2=geom.bounds[3]
    )
)


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
  • 11,011
  • 3
  • 34
  • 69
0

@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[1],
                lat_2=shape.bounds[3]
                ),
            always_xy=True
        )
David Ignjic
  • 1
  • 1