Área em KM do polígono de coordenadas

14

Eu tenho polígonos de coordenadas em (python shapely) que se parece com isso

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

Gostaria de calcular a área deste polígono em km ^ 2. Qual seria a melhor maneira de fazer isso em Python?

python polygon area shapely amaatouq
fonte
1
Você pode olhar para stackoverflow.com/questions/23697374/…
SIslam
Se você está recebendo o seguinte erro ao implementar uma das soluções acima, é porque lat1 e lat2 devem ser lat_1 e lat_2: pyproj.exceptions.CRSError: Projeção inválida: + proj = aea + lat1 = 37.843975868971484 + lat2 = 37.844325658890924 + type = crs: (Proj Erro interno: proj_create: Erro -21: cónica lat_1 = -lat_2)
Ramtin Kermani

Respostas:

20

Não estava prontamente aparente para mim como usar a resposta @sgillies, então aqui está uma versão mais detalhada:

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
fonte
1
O valor resultante não é exatamente o mesmo que o obtido em geojson.io . Por quê?
Astrojuanlu 06/07
16

Parece que suas coordenadas são longitude e latitude, sim? Use a shapely.ops.transformfunção Shapely para transformar o polígono em coordenadas projetadas de área igual e, em seguida, ocupe a área.

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
fonte
1
Você provavelmente deve indicar que partialnão é um built-in; pyprojterá que ser importado e, possivelmente instalado, etc.
kingledion
Eu notei que pyproj.Proj(proj='aea', lat1=lat1, lat2=lat2)resultou CRSError: Invalid projection: +proj=aea +lat1=5.0 +lat2=6.0 +type=crs. Alterar lat{1,2}em lat_{1,2}como implicado pela documentação proj4 fixa-lo: pyproj.Proj(proj='aea', lat1=lat1, lat2=lat2). Esta resposta é precisa ou deve ser atualizada?
Herbert
1
Eu precisava usar lat_1e em lat_2vez de lat1e lat2. Eu suspeito que isso se aplica post PROJ 6.0.0
oortCloud 26/04/19
3

As respostas acima parecem estar corretas, EXCETO que em algum momento recentemente, os parâmetros lat1 e lat2 no código pyproj foram renomeados com sublinhados: lat_1 e lat_2 (consulte /programming//a/55259718/1538758 ). Não tenho representante suficiente para comentar, por isso estou fazendo uma nova resposta (desculpe, desculpe)

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
fonte
1

Eu tropecei na "área", que parece mais simples de usar:

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

Por exemplo:

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))

... retorna:

área m2: 1082979.880942425

área km2: 1,082979880942425

Mike Honey
fonte