4

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'][0].area)
print(merged_df['geometry'][1].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?

enter image description here enter image description here

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?

2

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'][0].area)
print(merged_df['geometry'][1].area)

# OUTPUT: (in m2)
# 3098282528.696134
# 4325041659.243719
| improve this answer | |
  • 2
    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
  • 1
    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

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.