Fläche in KM vom Koordinatenpolygon

13

Ich habe Polygone aus Koordinaten (pythonförmig), die so aussehen

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

Ich möchte die Fläche dieses Polygons in km ^ 2 berechnen. Was wäre der beste Weg, um dies in Python zu tun?

python polygon area shapely amaatouq
quelle
1
Sie können sich stackoverflow.com/questions/23697374/…
SIslam
Wenn beim Implementieren einer der obigen Lösungen der folgende Fehler auftritt, sollten lat1 und lat2 lat_1 und lat_2 lauten: pyproj.exceptions.CRSError: Ungültige Projektion: + proj = aea + lat1 = 37.843975868971484 + lat2 = 37.844325658890924 + type = crs: (Interner Projektfehler: proj_create: Fehler -21: conic lat_1 = -lat_2)
Ramtin Kermani

Antworten:

19

Mir war nicht sofort klar, wie @sgillies answer zu verwenden ist. Hier ist eine ausführlichere Version:

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
quelle
1
Der resultierende Wert stimmt nicht genau mit dem in geojson.io ermittelten Wert überein . Warum?
Astrojuanlu
15

Es sieht so aus, als wären deine Koordinaten Längen- und Breitengrad, ja? Verwenden Sie die shapely.ops.transformFunktion von Shapely, um das Polygon in projizierte Koordinaten gleicher Fläche zu transformieren und dann die Fläche zu übernehmen.

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
quelle
1
Sie sollten wahrscheinlich angeben, dass partiales sich nicht um eine integrierte Funktion handelt. pyprojmüssen importiert und möglicherweise installiert werden, etc.
kingledion
Ich bemerkte, pyproj.Proj(proj='aea', lat1=lat1, lat2=lat2)führte zu CRSError: Invalid projection: +proj=aea +lat1=5.0 +lat2=6.0 +type=crs. Ändern lat{1,2}in lat_{1,2}wie implizierten proj4 Dokumentation es fest: pyproj.Proj(proj='aea', lat1=lat1, lat2=lat2). Ist diese Antwort korrekt oder sollte sie aktualisiert werden?
Herbert
1
Ich musste lat_1und lat_2anstelle von lat1und verwenden lat2. Ich vermute, dies gilt nach PROJ 6.0.0
oortCloud
2

Die obigen Antworten scheinen korrekt zu sein, AUSSER, dass die Parameter lat1 und lat2 im PyProj-Code kürzlich mit Unterstrichen umbenannt wurden: lat_1 und lat_2 (siehe /programming//a/55259718/1538758 ). Ich habe nicht genug Repräsentanten, um einen Kommentar abzugeben, also gebe ich eine neue Antwort (sorry, nicht sorry)

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

Ich bin über "Bereich" gestolpert, der einfacher zu bedienen zu sein scheint:

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

Beispielsweise:

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

... kehrt zurück:

Fläche m2: 1082979.880942425

Fläche km2: 1,082979880942425

Mike Honey
quelle