Área en KM desde Polígono de coordenadas

14

Tengo polígonos de coordenadas en (Python bien formados) que se ve así

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

Me gustaría calcular el área de este polígono en km ^ 2. ¿Cuál sería la mejor manera de hacer esto en Python?

python polygon area shapely amaatouq
fuente
1
Puede mirar stackoverflow.com/questions/23697374/…
SIslam
Si obtiene el siguiente error al implementar una de las soluciones anteriores, se debe a que lat1 y lat2 deberían ser lat_1 y lat_2: pyproj.exceptions.CRSError: proyección no válida: + proj = aea + lat1 = 37.843975868971484 + lat2 = 37.844325658890924 + tipo = crs: (Error interno del proyecto: proj_create: Error -21: conic lat_1 = -lat_2)
Ramtin Kermani

Respuestas:

20

No me resultó evidente cómo usar la respuesta @sgillies, así que aquí hay una versión más detallada:

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
fuente
1
El valor resultante no es exactamente el mismo que el obtenido en geojson.io . ¿Por qué?
astrojuanlu
16

Parece que tus coordenadas son longitud y latitud, ¿sí? Use la shapely.ops.transformfunción de Shapely para transformar el polígono en coordenadas de área igual proyectadas y luego tome el á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
fuente
1
Probablemente debería indicar que partialno está integrado; pyprojtendrá que ser importado y posiblemente instalado, etc
kingledion
Me di cuenta de que pyproj.Proj(proj='aea', lat1=lat1, lat2=lat2)resultó en CRSError: Invalid projection: +proj=aea +lat1=5.0 +lat2=6.0 +type=crs. Cambiando lat{1,2}a lat_{1,2}como se deduce de la documentación proj4 fija que: pyproj.Proj(proj='aea', lat1=lat1, lat2=lat2). ¿Es correcta esta respuesta o debería actualizarse?
Herbert
1
Necesitaba usar lat_1y en lat_2lugar de lat1y lat2. Sospecho que esto se aplica después de PROJ 6.0.0
oortCloud
3

Las respuestas anteriores parecen ser correctas, EXCEPTO que en algún momento recientemente, los parámetros lat1 y lat2 en el código pyproj fueron renombrados con guiones bajos: lat_1 y lat_2 (consulte /programming//a/55259718/1538758 ). No tengo suficiente representante para comentar, así que estoy haciendo una nueva respuesta (lo siento, no lo siento)

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

Me topé con "área" que parece más simple de usar:

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

Por ejemplo:

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

... devoluciones:

Área m2: 1082979.880942425

área km2: 1.082979880942425

Mike Honey
fuente