Search code examples
pythongisshapelypyproj

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

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)? enter image description here

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

Solution

  • 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