Mendapatkan area poligon menggunakan geopanda?

16

Diberikan geopandas GeoDataFrameserangkaian poligon, saya ingin mendapatkan luas dalam km persegi dari setiap fitur dalam daftar saya.

Ini adalah masalah yang cukup umum, dan solusi yang disarankan di masa lalu adalah menggunakan shapelydan pyprojsecara langsung (misalnya di sini dan di sini ).

Apakah ada cara untuk melakukan ini dengan murni geopandas?

Aleksey Bilogur
sumber

Jawaban:

17

Jika crs GeoDataFrame diketahui (EPSG: 4326 unit = degree, di sini), Anda tidak perlu Shapely, atau pyproj dalam skrip Anda karena GeoPandas menggunakannya).

import geopandas as gpd
test = gpd.read_file("test_wgs84.shp")
print test.crs
test.head(2)

masukkan deskripsi gambar di sini

Sekarang salin GeoDataFrame Anda dan ubah proyeksi ke sistem Cartesian (EPSG: 3857, unit = m seperti pada jawaban ResMar)

tost = test.copy()
tost= tost.to_crs({'init': 'epsg:3857'})
print tost.crs
tost.head(2)

masukkan deskripsi gambar di sini

Sekarang area dalam kilometer persegi

tost["area"] = tost['geometry'].area/ 10**6
tost.head(2)

masukkan deskripsi gambar di sini

Tetapi permukaan dalam proyeksi Mercator tidak benar, demikian juga dengan proyeksi lainnya dalam meter.

tost= tost.to_crs({'init': 'epsg:32633'})
tost["area"] = tost['geometry'].area/ 10**6
tost.head(2)

masukkan deskripsi gambar di sini

gen
sumber
Teks Anda benar epsg:3857, tetapi kode Anda benar epsg:3395, mana di antara keduanya yang benar?
Aleksey Bilogur
4
The .to_crsFungsi akan dilewatkan ke pyprojlagian. Contoh yang baik dari proyeksi area yang sama: proj4.org/projections/cea.html yang dapat diteruskan sebagai berikut:.to_crs({'proj':'cea'})
Swier
Untuk shapefile US Census Tracts setidaknya, saya dapat mengonfirmasi bahwa {'proj':'cea'}menghasilkan estimasi area terdekat.
Polor Beer
4

Saya percaya ya. Berikut ini harus bekerja:

gdf['geometry'].to_crs({'init': 'epsg:3395'})\
               .map(lambda p: p.area / 10**6)

Ini mengubah geometri menjadi proyeksi dengan luas yang sama, mengambil shapelyarea tersebut (dikembalikan dalam m ^ 2), dan memetakannya menjadi km ^ 2 (langkah terakhir ini opsional).

Aleksey Bilogur
sumber
Apakah ini benar?
Aleksey Bilogur
1
EPSG 3857 bukan area yang sama. en.wikipedia.org/wiki/Map_projection#Equal-area
alphabetasoup
Saya telah memodifikasi jawaban ini agar sesuai dengan epsg:3395CRS gen . Terima kasih.
Aleksey Bilogur