问题描述
我有GeoJSON
{
"type": "FeatureCollection",
"features": [
{
"type": "Feature",
"properties": {},
"geometry": {
"type": "Polygon",
"coordinates": [
[[13.65374516425911, 52.38533382814119], [13.65239769133293, 52.38675829106993], [13.64970274383571, 52.38675829106993], [13.64835527090953, 52.38533382814119], [13.64970274383571, 52.38390931824483], [13.65239769133293, 52.38390931824483], [13.65374516425911, 52.38533382814119]]
]
}
}
]
}
我想用Python计算其面积(87106.33m ^ 2).我该怎么办?
I would like to calculate its area (87106.33m^2) with Python. How do I do that?
# core modules
from functools import partial
# 3rd pary modules
from shapely.geometry import Polygon
from shapely.ops import transform
import pyproj
l = [[13.65374516425911, 52.38533382814119, 0.0], [13.65239769133293, 52.38675829106993, 0.0], [13.64970274383571, 52.38675829106993, 0.0], [13.64835527090953, 52.38533382814119, 0.0], [13.64970274383571, 52.38390931824483, 0.0], [13.65239769133293, 52.38390931824483, 0.0], [13.65374516425911, 52.38533382814119, 0.0]]
polygon = Polygon(l)
print(polygon.area)
proj = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'),
pyproj.Proj(init='epsg:3857'))
print(transform(proj, polygon).area)
它给出了1.1516745933889345e-05
和233827.03300877335
-第一个没有意义,但是我该如何解决第二个问题呢? (我不知道如何设置pyproj.Proj
init参数)
It gives 1.1516745933889345e-05
and 233827.03300877335
- that the first one doesn't make any sense was expected, but how do I fix the second one? (I have no idea how to set the pyproj.Proj
init parameter)
我想epsg:4326
是WGS84(源) ,但对于 epsg:3857 我不确定.
I guess epsg:4326
makes sense at it is WGS84 (source), but for epsg:3857 I'm uncertain.
以下内容更近了:
# core modules
from functools import partial
# 3rd pary modules
import pyproj
from shapely.geometry import Polygon
import shapely.ops as ops
l = [[13.65374516425911, 52.38533382814119, 0],
[13.65239769133293, 52.38675829106993, 0],
[13.64970274383571, 52.38675829106993, 0],
[13.64835527090953, 52.38533382814119, 0],
[13.64970274383571, 52.38390931824483, 0],
[13.65239769133293, 52.38390931824483, 0],
[13.65374516425911, 52.38533382814119, 0]]
polygon = Polygon(l)
print(polygon.area)
geom_area = ops.transform(
partial(
pyproj.transform,
pyproj.Proj(init='EPSG:4326'),
pyproj.Proj(
proj='aea',
lat1=polygon.bounds[1],
lat2=polygon.bounds[3])),
polygon)
print(geom_area.area)
它给出了87254.7m ^ 2-与geojson.io所说的仍然相差148m ^ 2.为什么会这样呢?
it gives 87254.7m^2 - that is still 148m^2 different from what geojson.io says. Why is that the case?
推荐答案
看起来geojson.io不在像您一样将球形坐标投影到平面上之后计算面积,而是使用一种特定的算法来计算面积直接从WGS84坐标在球体表面上绘制多边形的角度.如果要重新创建它,可以在此处找到源代码. >.
It looks like geojson.io is not calculating the area after projecting the spherical coordinates onto a plane like you are, but rather using a specific algorithm for calculating the area of a polygon on the surface of a sphere, directly from the WGS84 coordinates. If you want to recreate it you can find the source code here.
如果您乐意将坐标投影到平面系统以计算面积,因为它对于您的用例而言足够准确,那么您可以尝试使用此对德国的预测.例如:
If you are happy to carry on projecting the coordinates to a flat system to calculate the area, since it's good enough accuracy for your use case, then you might trying using this projection for Germany instead. E.g:
from osgeo import ogr
from osgeo import osr
source = osr.SpatialReference()
source.ImportFromEPSG(4326)
target = osr.SpatialReference()
target.ImportFromEPSG(5243)
transform = osr.CoordinateTransformation(source, target)
poly = ogr.CreateGeometryFromJson(str(geoJSON['features'][0]['geometry']))
poly.Transform(transform)
poly.GetArea()
返回87127.2534625642
这篇关于如何使用Python获取GeoJSON多边形的面积的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!