Polygonflächen mit Geopandas erhalten?

16

Angesichts einer geopandas GeoDataFrameReihe von Polygonen möchte ich die Fläche der einzelnen Elemente in meiner Liste in km² anzeigen.

Dies ist ein ziemlich häufiges Problem, und die in der Vergangenheit übliche Lösung bestand darin, shapelyund pyprojdirekt zu verwenden (z. B. hier und hier ).

Gibt es eine Möglichkeit, dies in rein zu tun geopandas?

Aleksey Bilogur
quelle

Antworten:

17

Wenn die crs des GeoDataFrame bekannt sind (EPSG: 4326 unit = degree, hier), benötigen Sie weder Shapely noch Pyproj in Ihrem Skript, da GeoPandas sie verwendet.

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

Bildbeschreibung hier eingeben

Kopieren Sie nun Ihren GeoDataFrame und ändern Sie die Projektion in ein kartesisches System (EPSG: 3857, Einheit = m wie in der Antwort von ResMar)

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

Bildbeschreibung hier eingeben

Jetzt ist die Fläche in Quadratkilometern

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

Bildbeschreibung hier eingeben

Die Flächen in der Mercator-Projektion sind jedoch nicht korrekt, also bei anderen Projektionen in Metern.

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

Bildbeschreibung hier eingeben

Gen
quelle
Ihr Text ist epsg:3857, aber Ihr Code ist epsg:3395, welcher der beiden ist korrekt?
Aleksey Bilogur
4
Die .to_crsFunktion wird pyprojtrotzdem weitergereicht. Ein gutes Beispiel für eine flächengleiche Projektion: proj4.org/projections/cea.html , die wie folgt übergeben werden kann:.to_crs({'proj':'cea'})
Swier
Zumindest für die Shapefiles der US-Volkszählungsdaten kann ich bestätigen, {'proj':'cea'}dass die nächstgelegenen Flächenschätzungen erstellt werden.
Polor Beer
4

Ich glaube ja Folgendes sollte funktionieren:

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

Dadurch wird die Geometrie in eine flächengleiche Projektion konvertiert, die Fläche abgerufen shapely(zurückgegeben in m ^ 2) und auf km ^ 2 abgebildet (dieser letzte Schritt ist optional).

Aleksey Bilogur
quelle
Ist das richtig?
Aleksey Bilogur
1
EPSG 3857 ist nicht gleich groß. en.wikipedia.org/wiki/Map_projection#Equal-area
alphabetasoup
Ich habe diese Antwort an das epsg:3395CRS des Gens angepasst. Vielen Dank.
Aleksey Bilogur