¿Conseguir áreas de polígonos usando geopandas?

16

Dado que geopandas GeoDataFramecontiene una serie de polígonos, me gustaría obtener el área en km² de cada entidad en mi lista.

Este es un problema bastante común, y la solución habitual sugerida en el pasado ha sido usar shapelyy pyprojdirectamente (por ejemplo, aquí y aquí ).

¿Hay alguna manera de hacer esto en puro geopandas?

Aleksey Bilogur
fuente

Respuestas:

17

Si se conocen los crs del GeoDataFrame (EPSG: 4326 unit = degree, aquí), no necesita Shapely, ni pyproj en su script porque GeoPandas los usa).

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

ingrese la descripción de la imagen aquí

Ahora copie su GeoDataFrame y cambie la proyección a un sistema cartesiano (EPSG: 3857, unit = m como en la respuesta de ResMar)

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

ingrese la descripción de la imagen aquí

Ahora el área en kilómetros cuadrados

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

ingrese la descripción de la imagen aquí

Pero las superficies en la proyección de Mercator no son correctas, así que con otras proyecciones en metros.

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

ingrese la descripción de la imagen aquí

gene
fuente
Su texto es epsg:3857, pero su código es epsg:3395, ¿cuál de los dos es correcto?
Aleksey Bilogur
44
La .to_crsfunción se pasa de pyprojtodos modos. Un buen ejemplo de una proyección de área igual: proj4.org/projections/cea.html que se puede aprobar de la siguiente manera:.to_crs({'proj':'cea'})
Swier
Por lo menos para los archivos de formas del Censo de los Estados Unidos, puedo confirmar que {'proj':'cea'}producen las estimaciones de área más cercanas.
Polor Beer
4

Yo creo que si. Lo siguiente debería funcionar:

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

Esto convierte la geometría en una proyección de área igual, recupera el shapelyárea (devuelta en m ^ 2) y la asigna a km ^ 2 (este último paso es opcional).

Aleksey Bilogur
fuente
¿Es esto correcto?
Aleksey Bilogur
1
EPSG 3857 no es área igual. en.wikipedia.org/wiki/Map_projection#Equal-area
alphabetasoup
He modificado esta respuesta para que se ajuste al epsg:3395CRS del gen . Gracias.
Aleksey Bilogur