Obtenir des zones polygonales à l'aide de géopandas?

16

Étant donné une geopandas GeoDataFramecontenant une série de polygones, je voudrais obtenir la superficie en km2 de chaque entité de ma liste.

C'est un problème assez courant, et la solution habituelle suggérée dans le passé a été d'utiliser shapelyet pyprojdirectement (par exemple ici et ici ).

Existe-t-il un moyen de le faire en pur geopandas?

Aleksey Bilogur
la source

Réponses:

17

Si les crs du GeoDataFrame sont connus (EPSG: 4326 unité = degré, ici), vous n'avez pas besoin de Shapely, ni pyproj dans votre script car GeoPandas les utilise).

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

entrez la description de l'image ici

Copiez maintenant votre GeoDataFrame et changez la projection en un système cartésien (EPSG: 3857, unité = m comme dans la réponse de ResMar)

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

entrez la description de l'image ici

Maintenant, la zone en kilomètres carrés

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

entrez la description de l'image ici

Mais les surfaces de la projection Mercator ne sont pas correctes, donc avec d'autres projections en mètres.

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

entrez la description de l'image ici

gène
la source
Votre texte est epsg:3857, mais votre code est epsg:3395, lequel des deux est correct?
Aleksey Bilogur
4
La .to_crsfonction est de pyprojtoute façon transmise à . Un bon exemple de projection à surface égale: proj4.org/projections/cea.html qui peut être passé comme suit:.to_crs({'proj':'cea'})
Swier
Pour les fichiers de formes des secteurs de recensement des États-Unis au moins, je peux confirmer qu'ils {'proj':'cea'}produisent les estimations de zone les plus proches.
Polor Beer du
4

Je crois que oui. Les éléments suivants devraient fonctionner:

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

Cela convertit la géométrie en une projection à surface égale, récupère la shapelyzone (renvoyée en m ^ 2) et la transforme en km ^ 2 (cette dernière étape est facultative).

Aleksey Bilogur
la source
Est-ce correct?
Aleksey Bilogur
1
EPSG 3857 n'est pas une surface égale. en.wikipedia.org/wiki/Map_projection#Equal-area
alphabetasoup
J'ai modifié cette réponse pour l'adapter au epsg:3395CRS du gène . Merci.
Aleksey Bilogur