【发布时间】：2021-06-24 15:45:38
【问题描述】：
我想找到根据纬度和经度创建的多边形的面积（以平方米为单位）。
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。
编辑：我交换了我的值，它必须是“经度”、“纬度”。感谢雪人2github
【问题讨论】：
-
你为什么不喜欢191？在任何情况下，您都不会告诉我们您正在使用哪个数据（“纬度/经度”坐标并不能告诉我们太多，因此图书馆不会理解它）。在输出上，你应该告诉它你想要米（你需要转换到“米坐标系”（通常用 x,y 表示）。在每一步你都应该写一个关于“基准”和“单位”的注释。这将有助于发现问题和丢失的转化。
-
@GiacomoCatenazzi 谢谢你的回复，我不想要
191，因为我的目标是获得 m^2 中的区域，并按照我帖子第一个链接中的代码，答案应该在 m^2 中，你能解释一下“你正在使用哪个数据”的意思吗，有没有一种方法可以告诉我正在使用纬度和经度？ -
纬度和经度只是调用两个角度的一种方式（例如使用 x,y），但它并没有真正定义精确的指向地球。您需要指定如何测量角度等（这是基准）。不同的国家（和不同的时间）使用不同的原点/偏移量（还有角度），还有地球的形状。所有地理程序都需要这样的数据（你在 pyproj 上给出了一个：WGS84（又名“GPS 坐标”）。数据是使用投影时最令人困惑的事情之一（我们都会犯很多错误）：所以要明确地写出来（并仔细检查）。
-
@GiacomoCatenazzi 好的，这对我来说很有意义，我认为这将是一项不需要太多领域知识的任务，但我错了