如何使用Python获取GeoJSON多边形的面积

Question

如何使用Python获取GeoJSON多边形的面积

12

我有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]]
        ]
      }
    }
  ]
}

这是一个关于IT技术的翻译任务，内容如下：

http://geojson.io 显示如下：

我想使用Python计算它的面积（87106.33平方米）。请问应该怎么做？

我的尝试

# 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的初始化参数）

我猜epsg:4326是有意义的，因为它是WGS84（来源），但对于epsg:3857我不确定。

更好的结果

以下结果更接近：

# 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.7平方米，与geojson.io给出的结果还有148平方米的差距。这是为什么呢？

- Martin Thoma

因为Web Mercator不能保留面积: http://www.geography.hunter.cuny.edu/~jochen/GTECH361/lectures/lecture04/concepts/Map%20coordinate%20systems/Map%20projections%20and%20distortion.htm。欢迎来到地图投影。你的问题更适合在[GIS.SE]上提问（你的*代码*没有问题，只是你选择的投影有问题），并且最好具体说明获取错误面积值的问题。这里是实际3857投影的链接：[https://epsg.io/3857](https://epsg.io/3857)，你提供的链接是其他投影。 - jpmc26

对于任何尝试使用附加的Python代码并出现错误的人，应该将lat1和lat2更改为lat_1和lat_2。 - Hannesh

2个回答

3

有一个适用于Python的Mapbox geojson-area端口：https://pypi.org/project/area/

from area import area

geoJSON = json.loads(...)

print(area(geoJSON['features'][0]['geometry']))

使用您的GeoJSON，它会返回87106.33，就像在http://geojson.io/上一样。

- Benedikt Köppel

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- fordy · Accepted Answer

看起来geojson.io在投影球面坐标到平面后没有计算面积，而是使用了一种特定的算法直接从WGS84坐标计算球面多边形的面积。如果您想重新创建它，可以在这里找到源代码。

如果您希望继续将坐标投影到平面系统中计算面积，在您的应用场景中足够准确，那么您可以尝试使用此德国投影方式。例如：

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。