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

Question

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

其中 http://geojson.io 显示为

我想用Python计算它的面积（87106.33m^2）。我该怎么做？

我尝试了什么

# 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.7m^2 - 这仍然与 geojson.io 所说的不同 148m^2 。为什么会这样？

Answer 1

看起来geojson.io并不是像你一样将球面坐标投影到平面上后计算面积，而是使用特定的算法直接从WGS84坐标计算球体表面上的多边形的面积。如果您想重新创建它，您可以在here找到源代码。

如果您愿意继续将坐标投影到平面系统来计算面积，因为它对于您的用例来说足够准确，那么您可以尝试使用 this 投影来代替德国。例如：

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

Answer 2

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

from area import area

geoJSON = json.loads(...)

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

使用您的 GeoJSON，它会返回

87106.33

，就像 http://geojson.io/。

Answer 3

如果使用GeoDjango，如here所述，您可以首先投影geojson并从中检索以平方米为单位的面积：

def get_acres(self): 
    """ 
    Returns the area in acres. 
    """ 
    # Convert our geographic polygons (in WGS84)
    # into a local projection - for New York here as an example (EPSG:32118)
    self.polygon.transform(32118)
    meters_sq = self.polygon.area
    
    acres = meters_sq * 0.000247105381  # meters^2 to acres

    return acres