Si le crs du GeoDataFrame est connu (EPSG:4326 unit=degree, ici), vous n'avez pas besoin de Shapely, ni de 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)
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)
Maintenant la superficie en kilomètres carrés
tost["area"] = tost['geometry'].area/ 10**6
tost.head(2)
Mais les surfaces dans la projection de 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)
Je crois que oui. Ce qui suit devrait 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 le shapely
zone (renvoyée en m^2), et la mappe à un km^2 (cette dernière étape est facultative).