问题描述
我想找到由纬度和经度创建的多边形的面积（以平方米为单位）。
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,] # 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解决方法
您可以使用 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,]
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