24

Given a GeoPandas's GeoDataFrame containing a series of polygons, I would like to get the area in km sq of each feature in my list.

This is a pretty common problem, and the usual suggested solution in the past has been to use shapely and pyproj directly (e.g. here and here).

Is there a way to do this in pure GeoPandas?

26

If the crs of the GeoDataFrame is known (EPSG:4326 unit=degree, here), you don't need Shapely, nor pyproj in your script because GeoPandas uses them).

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

enter image description here

Now copy your GeoDataFrame and change the projection to a Cartesian system (EPSG:3857, unit= m as in the answer of ResMar)

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

enter image description here

Now the area in square kilometers

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

enter image description here

But the surfaces in the Mercator projection are not correct, so with other projection in meters.

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

enter image description here

3
  • Your text is epsg:3857, but your code is epsg:3395, which of the two is correct? – Aleksey Bilogur Nov 20 '16 at 20:11
  • 6
    The .to_crs function gets passed to pyproj anyways. A good example of an equal area projection: proj4.org/projections/cea.html which can be passed as follows: .to_crs({'proj':'cea'}) – Swier Jan 24 '17 at 11:25
  • 1
    For the US Census Tracts shapefiles at least, I can confirm that {'proj':'cea'} produce the closest area estimations. – Polor Beer Dec 7 '17 at 20:12
8

I believe yes. The following ought to work:

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

This converts the geometry to an equal-area projection, fetches the shapely area (returned in m^2), and maps that to a km^2 (this last step is optional).

3
4

Yes, simply be sure to reproject your shape in Cylindrical equal-area format with {'proj':'cea'} that preserve area measure.

Then you can use .area method of your GeoDataFrame.

Your also need to divide by 1000000 because .area method give area in square meters.

import geopandas as gpd

gdf = gpd.read_file("YOUR_SHAPE_FILE.shp")
gdf = gdf['geometry'].to_crs({'proj':'cea'}) 

gdf.area / 10**6
0

Just a quick thought on the appropriate EPSG code for an equal-area estimation - 6933 may be a better "generic" solution (see https://epsg.io/6933 / https://www.mdpi.com/2220-9964/1/1/32)

No perfect solution for obvious reasons, but 6933 does seem to nicely merge the benefits of a cea and Lambert equal area.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.