python - Shapely/Pyproj 查找根据纬度和经度创建的多边形的面积(以 m^2 为单位)

Question

我想找到根据纬度和经度创建的多边形的面积(以平方米为单位)。

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

值191.56938242734225不是我想要的，所以我做了some search online，发现我必须转换并使用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

但是我得到nan，我在这里做错了什么？根据上面链接上的 a comment，对于我确实可以正常工作的某些数据(对于不同的纬度经度列表)，它与图像中显示的区域不匹配。

如何获取此图像中所示的区域(来自 geojson.io)？

编辑(尝试解决方案后)

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)

我仍然得到nan。

编辑:我交换了我的值，它必须是“经度”，“纬度”。感谢 github 的 Snowman2

Answer 1

您可以使用 pyproj 获取测地面积: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

另请参阅:https://pyproj4.github.io/pyproj/stable/api/geod.html#pyproj.Geod.geometry_area_perimeter

您可能需要shapely.ops.orient