Площадь в км от полигона координат

14

У меня есть полигоны из координат в (Python shapely), который выглядит следующим образом

POLYGON ((24.8085317 46.8512821, 24.7986952 46.8574619, 24.8088238 46.8664741, 24.8155239 46.8576335, 24.8085317 46.8512821))

Я хотел бы рассчитать площадь этого многоугольника в км ^ 2. Что было бы лучшим способом сделать это в Python?

amaatouq
источник
1
Вы можете посмотреть на stackoverflow.com/questions/23697374/...
SIslam
Если вы получаете следующую ошибку при реализации одного из вышеуказанных решений, это потому, что lat1 и lat2 должны быть lat_1 и lat_2: pyproj.exceptions.CRSError: Неверная проекция: + proj = aea + lat1 = 37.843975868971484 + lat2 = 37.844325658890924 + type = crs: (Внутренняя ошибка Proj: proj_create: Ошибка -21: conic lat_1 = -lat_2)
Рамтин Кермани

Ответы:

20

Мне было не совсем понятно, как использовать ответ @sgillies, поэтому вот более подробная версия:

import pyproj    
import shapely
import shapely.ops as ops
from shapely.geometry.polygon import Polygon
from functools import partial


geom = Polygon([(0, 0), (0, 10), (10, 10), (10, 0), (0, 0)])
geom_area = ops.transform(
    partial(
        pyproj.transform,
        pyproj.Proj(init='EPSG:4326'),
        pyproj.Proj(
            proj='aea',
            lat1=geom.bounds[1],
            lat2=geom.bounds[3])),
    geom)

# Print the area in m^2
print geom_area.area 
jczaplew
источник
1
Полученное значение не совсем совпадает с полученным в geojson.io . Почему?
astrojuanlu
16

Похоже, ваши координаты - это долгота и широта, да? Используйте shapely.ops.transformфункцию Shapely, чтобы преобразовать многоугольник в проекции равных координат области, а затем взять область.

python
import pyproj
from functools import partial

geom_aea = transform(
partial(
    pyproj.transform,
    pyproj.Proj(init='EPSG:4326'),
    pyproj.Proj(
        proj='aea',
        lat1=geom.bounds[1],
        lat2=geom.bounds[3])),
geom)

print(geom_aea.area)
# Output in m^2: 1083461.9234313113 
sgillies
источник
1
Вы должны вероятно указать, что partialэто не встроенный; pyprojпридется импортировать и, возможно, устанавливать и т. д.
kingledion
Я заметил, pyproj.Proj(proj='aea', lat1=lat1, lat2=lat2)что в результате CRSError: Invalid projection: +proj=aea +lat1=5.0 +lat2=6.0 +type=crs. Изменение lat{1,2}в lat_{1,2}качестве подразумеваемых proj4 документации установил ее: pyproj.Proj(proj='aea', lat1=lat1, lat2=lat2). Является ли этот ответ точным или должен быть обновлен?
Герберт
1
Мне нужно было использовать lat_1и lat_2вместо lat1и lat2. Я подозреваю, что это применимо пост PROJ 6.0.0
oortCloud
3

Приведенные выше ответы кажутся правильными, за исключением того, что в последнее время параметры lat1 и lat2 в коде pyproj были переименованы с подчеркиванием: lat_1 и lat_2 (см. Https://stackoverflow.com/a/55259718/1538758 ). У меня недостаточно представителей, чтобы комментировать, поэтому я делаю новый ответ (извините, не извините)

import pyproj    
import shapely
import shapely.ops as ops
from shapely.geometry.polygon import Polygon
from functools import partial


geom = Polygon([(0, 0), (0, 10), (10, 10), (10, 0), (0, 0)])
geom_area = ops.transform(
    partial(
        pyproj.transform,
        pyproj.Proj(init='EPSG:4326'),
        pyproj.Proj(
            proj='aea',
            lat_1=geom.bounds[1],
            lat_2=geom.bounds[3])),
    geom)

# Print the area in m^2
print geom_area.area 
MartinThurn
источник
1

Я наткнулся на «область», которая кажется проще в использовании:

https://pypi.org/project/area/

Например:

from area import area

obj = {'type':'Polygon','coordinates':[[[24.8085317,46.8512821], [24.7986952,46.8574619], [24.8088238,46.8664741], [24.8155239,46.8576335], [24.8085317,46.8512821]]]}

area_m2 = area(obj)

area_km2 = area_m2 / 1e+6
print ('area m2:' + str(area_m2))
print ('area km2:' + str(area_km2))

... возвращает:

площадь м2: 1082979.880942425

площадь км2: 1,082979880942425

Майк Хани
источник