# What is the area of my shapes?

I created overlapping geometries and I would like to know the area of the intersection.

``````import geopandas as gpd
from shapely.geometry import Polygon

polys1 = gpd.GeoSeries([Polygon([(0,0), (0,1), (1,1), (1,0)]),
Polygon([(89,89), (89,90), (90,90), (90,89)]),

])

polys2 = gpd.GeoSeries([Polygon([(0.5,0.5), (0.5,1.5), (1.5,1.5), (1.5,0.5)]),
Polygon([(88.5,88.5), (88.5,89.5), (89.5,89.5), (89.5,88.5)]),
])

df1 = gpd.GeoDataFrame({'geometry': polys1, 'seed_value':[10,10]})
df2 = gpd.GeoDataFrame({'geometry': polys2})

merged_df = gpd.overlay(df2, df1, how='intersection')

print(merged_df['geometry'].area)
print(merged_df['geometry'].area)
``````

One intersection is located at the equator (left) and the other at a pole (right). The code returns an area for both intersections of 0.25.

Shouldn't the area be different because of the different length of arcminutes at the equator and at the poles? How should I interpret the result and how can I compute the 'real' area?

If I would like to pass on values from one geometry to another geometry weighted by the overlapping area, what values do I have to compute to take into account the different lengths of arc-minutes at the poles and at the equator?

All geometric operations that GeoPandas performs depends on Shapely package. Shapely performs all operations in the x-y plane. That means Shapely doesn't care if a shape is located at equator or at a pole, or it has geographic or projected coordinate system.

0.25 means 0.25 squared degree (0.25 deg2) which is meaningless in terms of area calculation. If they were in a projected coordinate system, it could be a meaningful value (0.25 m2/km2/inch2/etc.)

For further information, please review: Coordinate Systems - Geometric Objects in Shapely

How to compute the 'real' area?

Actually, there has never been the real area in reality. It depends on coordinate system that you use. If you want to use a projected coordinate system (PCS), you should know that every PCS has distortions. Even if you use equal-area PCS, area depends on Geographic CS which PCS is based-on. It also depends on assuming if the earth is a sphere or ellipsoid.

But you can follow steps to get the projected area:

• First, you should set a crs for `df1` and `df2`:

``````# I assume both are in WGS 84
df1.crs = "EPSG:4326"
df2.crs = "EPSG:4326"
``````
• And change the projection to a Cartesian system, for example `EPSG:3857`(Web Mercator projection). As @Vince stated, Web Mercator is useless for global area calculation. I just use it as an example.

``````df1 = df1.to_crs(3857) # Web Mercator (unit: meter)
df2 = df2.to_crs(3857)
``````
• Then apply intersection:

``````merged_df = gpd.overlay(df2, df1, how='intersection')
``````

But notice that, in your case, you get an error. `IllegalArgumentException: RobustDeterminant encountered non-finite numbers`. Because in Web Mercator projection, coordinate value at a pole is infinite. Therefore, you should select a suitable projection for both equator and pole areas.

Results for different `polys1` and `polys2`:

``````import geopandas as gpd
from shapely.geometry import Polygon

polys1 = gpd.GeoSeries([Polygon([(0,0), (0,1), (1,1), (1,0)]),
Polygon([(44,44), (44,45), (45,45), (45,44)])])

polys2 = gpd.GeoSeries([Polygon([(0.5,0.5), (0.5,1.5), (1.5,1.5), (1.5,0.5)]),
Polygon([(43.5,43.5), (43.5,44.5), (44.5,44.5), (44.5,43.5)])])

df1 = gpd.GeoDataFrame({'geometry': polys1, 'seed_value':[10,10]})
df2 = gpd.GeoDataFrame({'geometry': polys2})

df1.crs = "EPSG:4326" # WGS 84 (unit: degree)
df2.crs = "EPSG:4326"

df1 = df1.to_crs(3857) # Web Mercator (unit: meter)
df2 = df2.to_crs(3857)

merged_df = gpd.overlay(df2, df1, how='intersection')

print(merged_df['geometry'].area)
print(merged_df['geometry'].area)

# OUTPUT: (in m2)
# 3098282528.696134
# 4325041659.243719
``````
• Careful! Web Mercator is just as useless for global area calculation as Cartesian degrees, possibly more so, since the nominal units are square meters but the values are nearly infinitely wrong. – Vince Feb 25 at 12:49
• @Vince Thank you for your warning. I'm just giving an example. – Kadir Şahbaz Feb 25 at 12:52
• @Vince, please feel free to edit the answer if there is an incorrect information. I answered as far as I know. – Kadir Şahbaz Feb 25 at 15:26
• Thank you vey very much for your detailed answer. So if I would like to distribute values from one GDF to another GDF based on overlying area, does it make sense at all to use GDF['geometry'].area in degrees or will that return wrong results? – Stücke Feb 25 at 15:34
• If you use a projected CS and it's sutiable for area you work, probably yes. You can use `gdf.area` instead of `gdf["geometry"].area`. They are the same. – Kadir Şahbaz Feb 25 at 15:49