Pozyskiwanie obszarów wielokątów za pomocą geopandaz?

16

Biorąc pod uwagę geopandas GeoDataFrameserię wielokątów, chciałbym uzyskać powierzchnię w km2 każdej funkcji na mojej liście.

Jest to dość powszechny problem, a zwykle proponowane rozwiązanie w przeszłości był w użyciu shapelyi pyprojbezpośrednio (np tutaj i tutaj ).

Czy istnieje sposób, aby to zrobić w czystej postaci geopandas?

Aleksey Bilogur
źródło

Odpowiedzi:

17

Jeśli znany jest crs GeoDataFrame (EPSG: 4326 jednostka = stopień, tutaj), nie potrzebujesz Shapely ani pyproj w swoim skrypcie, ponieważ używa go GeoPandas).

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

wprowadź opis zdjęcia tutaj

Teraz skopiuj GeoDataFrame i zmień rzut na system kartezjański (EPSG: 3857, jednostka = m jak w odpowiedzi ResMar)

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

wprowadź opis zdjęcia tutaj

Teraz obszar w kilometrach kwadratowych

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

wprowadź opis zdjęcia tutaj

Ale powierzchnie w rzucie Mercatora są niepoprawne, więc w przypadku innych rzutów w metrach.

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

wprowadź opis zdjęcia tutaj

gen
źródło
Twój tekst jest epsg:3857, ale kod brzmi epsg:3395, który z tych dwóch jest poprawny?
Aleksey Bilogur
4
Tak czy inaczej .to_crsfunkcja zostaje przekazana pyproj. Dobry przykład rzutu równego obszaru: proj4.org/projections/cea.html, który można przekazać w następujący sposób:.to_crs({'proj':'cea'})
Świer
Przynajmniej dla plików kształtowych Traktu ze Spisu Powszechnego, mogę potwierdzić, że {'proj':'cea'}dają one najbliższe oszacowania powierzchni.
Polor Beer,
4

Wierzę tak. Następujące powinny działać:

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

Konwertuje geometrię na rzut o jednakowym obszarze, pobiera shapelyobszar (zwrócony wm ^ 2) i odwzorowuje go na km ^ 2 (ten ostatni krok jest opcjonalny).

Aleksey Bilogur
źródło
Czy to jest poprawne?
Aleksey Bilogur
1
EPSG 3857 nie jest równym obszarem. en.wikipedia.org/wiki/Map_projection#Equal-area
alphabetasoup
Zmodyfikowałem tę odpowiedź, aby dopasować epsg:3395CRS genu . Dzięki.
Aleksey Bilogur