Area dalam KM dari Poligon koordinat

14

Saya memiliki poligon dari koordinat di (python rupawan) yang terlihat seperti ini

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

Saya ingin menghitung luas poligon ini dalam km ^ 2. Apa cara terbaik untuk melakukan ini dengan Python?

python polygon area shapely amaatouq
sumber
1
Anda dapat melihat stackoverflow.com/questions/23697374/…
SIslam
Jika Anda mendapatkan kesalahan berikut menerapkan salah satu solusi di atas, adalah karena lat1 dan lat2 harus lat_1 dan lat_2: pyproj.exceptions.CRSError: Proyeksi tidak valid: + proj = aea + lat1 = 37.843975868971484 + lat2 = 37.844325658890924 + tipe = crs: (Kesalahan Proj Internal: proj_create: Kesalahan -21: conic lat_1 = -lat_2)
Ramtin Kermani

Jawaban:

20

Bagi saya tidak jelas bagaimana cara menggunakan jawaban @sgillies, jadi di sini adalah versi yang lebih verbose:

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
sumber
1
Nilai yang dihasilkan tidak persis sama dengan yang diperoleh di geojson.io . Mengapa?
astrojuanlu
16

Sepertinya koordinat Anda garis bujur dan lintang, ya? Gunakan shapely.ops.transformfungsi Shapely untuk mengubah poligon untuk memproyeksikan koordinat area yang sama dan kemudian mengambil area tersebut.

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
sumber
1
Anda mungkin harus menunjukkan bahwa partialitu bukan bawaan; pyprojharus diimpor dan mungkin dipasang, dll.
kingledion
Saya perhatikan pyproj.Proj(proj='aea', lat1=lat1, lat2=lat2)menghasilkan CRSError: Invalid projection: +proj=aea +lat1=5.0 +lat2=6.0 +type=crs. Mengubah lat{1,2}ke lat_{1,2}seperti yang tersirat oleh dokumentasi PROJ4 tetap itu: pyproj.Proj(proj='aea', lat1=lat1, lat2=lat2). Apakah Jawaban ini akurat, atau haruskah diperbarui?
Herbert
1
Saya perlu menggunakan lat_1dan lat_2bukannya lat1dan lat2. Saya menduga ini berlaku post PROJ 6.0.0
oortCloud
3

Jawaban di atas tampaknya benar, KECUALI bahwa pada beberapa titik baru-baru ini, parameter lat1 dan lat2 dalam kode pyproj diganti nama dengan garis bawah: lat_1 dan lat_2 (lihat /programming//a/55259718/1538758 ). Saya tidak punya cukup perwakilan untuk berkomentar, jadi saya membuat jawaban baru (maaf bukan maaf)

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
Martin
sumber
1

Saya menemukan "area" yang tampaknya lebih mudah digunakan:

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

Sebagai contoh:

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

... pengembalian:

area m2: 1082979.880942425

area km2: 1.082979880942425

Mike Honey
sumber