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?

— 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

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