Aire en KM du polygone de coordonnées

14

J'ai des polygones de coordonnées en (python galbé) qui ressemble à ceci

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

Je voudrais calculer l'aire de ce polygone en km ^ 2. Quelle serait la meilleure façon de le faire en Python?

amaatouq
la source
1
Vous pouvez consulter stackoverflow.com/questions/23697374/…
SIslam
Si vous obtenez l'erreur suivante lors de l'implémentation de l'une des solutions ci-dessus, c'est parce que lat1 et lat2 doivent être lat_1 et lat_2: pyproj.exceptions.CRSError: Projection non valide: + proj = aea + lat1 = 37.843975868971484 + lat2 = 37.844325658890924 + type = crs: (Erreur Proj interne: proj_create: Erreur -21: conic lat_1 = -lat_2)
Ramtin Kermani

Réponses:

20

Il n'était pas évident pour moi comment utiliser la réponse @sgillies, alors voici une version plus détaillée:

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
la source
1
La valeur résultante n'est pas exactement la même que celle obtenue dans geojson.io . Pourquoi?
astrojuanlu
16

Il semble que vos coordonnées soient la longitude et la latitude, oui? Utilisez la shapely.ops.transformfonction de Shapely pour transformer le polygone en coordonnées de zone égales projetées , puis prenez la zone.

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
la source
1
Vous devriez probablement indiquer qu'il partialne s'agit pas d'un système intégré; pyprojdevra être importé et éventuellement installé, etc.
kingledion
J'ai remarqué un pyproj.Proj(proj='aea', lat1=lat1, lat2=lat2)résultat CRSError: Invalid projection: +proj=aea +lat1=5.0 +lat2=6.0 +type=crs. La modification lat{1,2}en lat_{1,2}tant que sous - entendus par la documentation PROJ4 fixe il: pyproj.Proj(proj='aea', lat1=lat1, lat2=lat2). Cette réponse est-elle exacte ou doit-elle être mise à jour?
Herbert
1
J'avais besoin d'utiliser lat_1et lat_2au lieu de lat1et lat2. Je soupçonne que cela s'applique après PROJ 6.0.0
oortCloud
3

Les réponses ci-dessus semblent être correctes, SAUF qu'à un moment donné récemment, les paramètres lat1 et lat2 dans le code pyproj ont été renommés avec des traits de soulignement: lat_1 et lat_2 (voir /programming//a/55259718/1538758 ). Je n'ai pas assez de représentant pour commenter, donc je fais une nouvelle réponse (désolé pas désolé)

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
la source
1

Je suis tombé sur "zone" qui semble plus simple à utiliser:

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

Par exemple:

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

... Retour:

surface m2: 1082979.880942425

superficie en km2: 1.082979880942425

Mike Honey
la source