Как получить площадь полигона GeoJSON с помощью Python

У меня есть 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 отображается как

введите описание изображения здесь

Я хотел бы рассчитать его площадь (87106,33 м ^ 2) с помощью Python. Как я могу это сделать?

Что я пробовал

# 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 м ^ 2 - это все еще 148 м ^ 2 отличается от того, что говорит geojson.io. Почему это так?


person Martin Thoma    schedule 27.07.2018    source источник
comment
Поскольку Web Mercator не сохраняет площадь: geography.hunter.cuny.edu/~jochen/GTECH361/lectures/lecture04/. Добро пожаловать в картографические проекции. Ваш вопрос лучше подходит для географических информационных систем (в вашем коде нет ничего плохого, просто ваш выбор проекции. ) и лучше указать конкретное значение для площади. И вот ссылка на реальную проекцию 3857; тот, на который вы ссылаетесь, - это что-то другое.   -  person jpmc26    schedule 27.07.2018


Ответы (1)


Похоже, что 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

person fordy    schedule 27.07.2018
comment
Хм... Есть ли простой ответ, что такое правильная проекция? Каковы основные различия между 4326 и 5243? - person Martin Thoma; 27.07.2018
comment
@MartinThoma Нет такой вещи, как простой ответ на правильную проекцию. Все проекции имеют разные свойства. Это очень сильно зависит от вашего конкретного варианта использования. В частности, многие прогнозы дают точные значения только в определенных регионах мира. - person jpmc26; 27.07.2018
comment
Это супер интересно! Не могли бы вы назвать два разумных прогноза для Берлина и их преимущества по сравнению друг с другом? - person Martin Thoma; 27.07.2018