Obtendo áreas de polígonos usando geopandas?

16

Dado que geopandas GeoDataFramecontém uma série de polígonos, gostaria de obter a área em km² de cada elemento da minha lista.

Este é um problema bastante comum, e a solução sugerida no passado foi usar shapelye pyprojdiretamente (por exemplo, aqui e aqui ).

Existe uma maneira de fazer isso de maneira pura geopandas?

Aleksey Bilogur
fonte

Respostas:

17

Se os CRs do GeoDataFrame são conhecidos (EPSG: 4326 unit = degree, aqui), você não precisa do Shapely, nem do pyproj em seu script, porque o GeoPandas os usa).

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

insira a descrição da imagem aqui

Agora copie seu GeoDataFrame e altere a projeção para um sistema cartesiano (EPSG: 3857, unit = m, como na resposta de ResMar)

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

insira a descrição da imagem aqui

Agora a área em quilômetros quadrados

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

insira a descrição da imagem aqui

Mas as superfícies na projeção Mercator não estão corretas, assim como outras projeções em metros.

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

insira a descrição da imagem aqui

gene
fonte
Seu texto é epsg:3857, mas seu código é epsg:3395, qual dos dois está correto?
Aleksey Bilogur
4
A .to_crsfunção é passada para pyprojqualquer maneira. Um bom exemplo de uma projeção de área igual: proj4.org/projections/cea.html, que pode ser transmitida da seguinte forma:.to_crs({'proj':'cea'})
Swier
Para os shapefiles dos EUA, pelo menos, posso confirmar que {'proj':'cea'}produzem as estimativas de área mais próximas.
Polor Beer
4

Eu acredito que sim O seguinte deve funcionar:

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

Isso converte a geometria em uma projeção de área igual, busca a shapelyárea (retornada em m ^ 2) e mapeia para um km ^ 2 (essa última etapa é opcional).

Aleksey Bilogur
fonte
Isso está correto?
274168 Aleksey Bilogur #
1
O EPSG 3857 não é uma área igual. pt.wikipedia.org/wiki/Map_projection#Equal-area
alphabetasoup
Modifiquei esta resposta para ajustar-se à epsg:3395CRS do gene . Obrigado.
Aleksey Bilogur