Wednesday, 7 August 2019

python - Getting polygon areas using geopandas?


Given a geopandas 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?



Answer



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


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...