Obszar w KM od wielokąta o współrzędnych

14

Mam wielokąty ze współrzędnych w (kształtny python), które wyglądają tak

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

Chciałbym obliczyć powierzchnię tego wielokąta w km ^ 2. Jaki byłby najlepszy sposób, aby to zrobić w Pythonie?

python polygon area shapely amaatouq
źródło
1
Możesz spojrzeć na stackoverflow.com/questions/23697374/…
SIslam
Jeśli pojawia się następujący błąd przy wdrażaniu jednego z powyższych rozwiązań, to dlatego, że lat1 i lat2 powinny być lat_1 i lat_2: pyproj.exceptions.CRSError: Niepoprawna projekcja: + proj = aea + lat1 = 37.843975868971484 + lat2 = 37.844325658890924 + typ = crs: (Internal Proj Error: proj_create: Error -21: conic lat_1 = -lat_2)
Ramtin Kermani

Odpowiedzi:

20

Nie było dla mnie oczywiste, jak korzystać z odpowiedzi @sgillies, więc oto bardziej szczegółowa wersja:

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
źródło
1
Wynikowa wartość nie jest dokładnie taka sama jak wartość uzyskana w geojson.io . Dlaczego?
astrojuanlu
16

Wygląda na to, że twoje współrzędne to długość i szerokość geograficzna, tak? Użyj shapely.ops.transformfunkcji Shapely, aby przekształcić wielokąt w rzutowane współrzędne równego obszaru, a następnie obierz obszar.

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
źródło
1
Prawdopodobnie powinieneś wskazać, że partialnie jest to wbudowane; pyprojbędą musiały zostać zaimportowane i prawdopodobnie zainstalowane itp.
kingledion
Zauważyłem, pyproj.Proj(proj='aea', lat1=lat1, lat2=lat2)że w rezultacie CRSError: Invalid projection: +proj=aea +lat1=5.0 +lat2=6.0 +type=crs. Zmiana lat{1,2}w lat_{1,2}jak sugeruje dokumentacja proj4 ustalone to: pyproj.Proj(proj='aea', lat1=lat1, lat2=lat2). Czy ta odpowiedź jest dokładna, czy powinna zostać zaktualizowana?
Herbert
1
Musiałem użyć lat_1i lat_2zamiast lat1i lat2. Podejrzewam, że dotyczy to postu PROJ 6.0.0
oortCloud,
3

Powyższe odpowiedzi wydają się poprawne, Z WYJĄTKIEM, że w pewnym momencie parametry lat1 i lat2 w kodzie pyproj zostały przemianowane na podkreślenia: lat_1 i lat_2 (patrz /programming//a/55259718/1538758 ). Nie mam wystarczającej liczby przedstawicieli, aby móc komentować, więc robię nową odpowiedź (przepraszam, nie przepraszam)

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
źródło
1

Natknąłem się na „obszar”, który wydaje się prostszy w użyciu:

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

Na przykład:

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

... zwroty:

powierzchnia m2: 1082979.880942425

powierzchnia km2: 1,082979880942425

Mike Honey
źródło