Shapely / Pyproj find area (in m^2) of a polygon created from latitude and longitude

Question

I want to find the area in square metre for a polygon created from latitude and longitude.

import shapely.ops as ops
from shapely.geometry import Polygon
from pyproj import transform, Proj
from functools import partial

value = [
    [83.3203125, 58.26328705248601],
    [98.7890625, 58.81374171570782],
    [105.1171875, 60.930432202923335],
    [104.0625, 65.6582745198266],
    [97.734375, 67.47492238478702],
    [87.890625, 67.06743335108298],
    [79.8046875, 65.36683689226321],
    [79.1015625, 62.431074232920906],
    [83.3203125, 58.26328705248601],
] # from geojson

polygon = Polygon(value)

print(polygon.area)  # 191.56938242734225

The value 191.56938242734225 is not what I want, so I did some search online and found out I had to transform and use pyproj.

area = ops.transform(
    partial(transform, Proj("EPSG:4326"), Proj(proj="aea", lat_1=polygon.bounds[1], lon_1=polygon.bounds[3])), polygon
)

print(area.area)  # nan

But I get nan, what am I doing wrong here? As per a comment on the link above, for some of the data where I do get this working (for a different latitude longitude list) it does not match with the area shown in the image.

How can I get the area as seen in this image (from geojson.io)?

Edit (after trying the solution)

from pyproj import Geod
from shapely.geometry import Polygon
from shapely.ops import orient

value = [
    (49.16848657161892, -122.4927565020954),
    (49.17145542056141, -122.44026815784645),
    (49.168944059578955, -122.42223809616848),
    (49.18182618668267, -122.42070673842014),
    (49.18901155779426, -122.40938116907655),
    (49.19333350498177, -122.41092556489616),
    (49.19754283081477, -122.4063781772052),
    (49.20350953892491, -122.38609910402336),
    (49.21439468669444, -122.36070237276812),
    (49.222360646955416, -122.35469638902534),
    (49.22411415726248, -122.35687285265936),
    (49.22498357252119, -122.35828854882728),
    (49.226469956724024, -122.35790244987238),
    (49.22832086322251, -122.35635805405282),
    (49.22936323360339, -122.35609949522244),
    (49.22913888972668, -122.35764389104204),
    (49.22773671741685, -122.35871638813896),
    (49.226739672517496, -122.36025878439766),
    (49.22010672437924, -122.38006780577793),
    (49.21996732086788, -122.38057257126464),
    (49.21999536911368, -122.39299208764706),
    (49.21992686126862, -122.4036494466543),
    (49.21891709518424, -122.40613785636424),
    (49.21888904632648, -122.4073390531128),
    (49.219562214519144, -122.40965564684215),
    (49.2199548917306, -122.41051364451972),
    (49.22074023679359, -122.4168628273335),
    (49.22121704735026, -122.421925013631),
    (49.22035066674655, -122.42482694635817),
    (49.21995799267972, -122.43289212452709),
    (49.22046285876402, -122.5120853101642),
    (49.219966402618745, -122.55650262699096),
    (49.255967240082846, -122.55753222420404),
    (49.256303572742965, -122.59614211969338),
    (49.27154824650424, -122.59648531876442),
    (49.3008416630168, -122.59480957171624),
    (49.300169602582415, -122.6655085803457),
    (49.29456874254022, -122.66036059428048),
    (49.28941538916121, -122.6675677747718),
    (49.288389695754155, -122.68427714208296),
    (49.27998615908695, -122.6956027114265),
    (49.25981182761499, -122.7177390515071),
    (49.25420638251272, -122.72099944268174),
    (49.241760018963866, -122.73953219251663),
    (49.227520721395656, -122.77317064879188),
    (49.22633007199472, -122.7634995116452),
    (49.21729887187899, -122.75397573742444),
    (49.20512626699452, -122.70094375077176),
    (49.19625956267232, -122.6707422325223),
    (49.19968297235687, -122.65383967827474),
    (49.210563761916895, -122.62311572400029),
    (49.21028322521899, -122.60552677161068),
    (49.19894821132008, -122.58853841759532),
    (49.1876984184996, -122.58114408968196),
    (49.16741121321411, -122.52297464048324),
]

polygon = Polygon(value)
geod = Geod(ellps="WGS84")
poly_area, poly_perimeter = geod.geometry_area_perimeter(orient(polygon))

print(poly_area)

I still get nan for this.

EDIT : I had my values swapped, it must be "longitude", "latitude". Thanks to snowman2 at github

Answer 1

You could use pyproj to get the geodesic area: https://pyproj4.github.io/pyproj/stable/examples.html#geodesic-area

from pyproj import Geod
from shapely.geometry import Polygon


value = [
    [83.3203125, 58.26328705248601],
    [98.7890625, 58.81374171570782],
    [105.1171875, 60.930432202923335],
    [104.0625, 65.6582745198266],
    [97.734375, 67.47492238478702],
    [87.890625, 67.06743335108298],
    [79.8046875, 65.36683689226321],
    [79.1015625, 62.431074232920906],
    [83.3203125, 58.26328705248601],
]

polygon = Polygon(value)
geod = Geod(ellps="WGS84")
poly_area, poly_perimeter = geod.geometry_area_perimeter(polygon)

print(poly_area)  # 1073367471345.5327

Also see: https://pyproj4.github.io/pyproj/stable/api/geod.html#pyproj.Geod.geometry_area_perimeter

You may need shapely.ops.orient