# Calculate intersection area in square meters of two polygons that use EPSG:4326 in Python

I have EPSG:4326 data for polygons and would like to calculate the intersection of these polygons in square meters.

Is it possible to do this without any transformations from EPSG:4326 to some other CRS like EPSG:3857?

``````def get_intersection(polygon1, polygon2):
area = 0.0
if polygon1.intersects(polygon2):
try:
area = polygon1.intersection(polygon2).area
print("area per units?? ", area)
except:
area = area
return area

if __name__ == "__main__":

poly1 = [(59.116666668946074, 5.133333344256348), (59.133333335610615, 5.1333333442660845), (59.133333335572246, 5.173333343890755), (59.13333333553494, 5.2133333435264975), (59.13333333549889, 5.253333343172972), (59.13333333546384, 5.293333342829971), (59.13333333542988, 5.333333342497181), (59.11666666876498, 5.333333342488369), (59.116666668799034, 5.293333342820962), (59.11666666883406, 5.253333343163808), (59.11666666887028, 5.213333343517124), (59.11666666890758, 5.1733333438812155)]

poly2 = [(59.127238679774464, 5.279128591100007), (59.127459966404096, 5.279065293749611), (59.127492503591576, 5.279495778853107), (59.12771379045561, 5.279432483234871), (59.127843926177974, 5.281154444053249), (59.12762263817868, 5.2812177288924875), (59.127655168284726, 5.281648218808487), (59.12699130313883, 5.281838059468042), (59.12686117762246, 5.280116137915042), (59.127082464918686, 5.280052847878851), (59.12704992971515, 5.2796223672021), (59.127271216678096, 5.279559073508844)]

polygon1 = Polygon(poly1)
polygon2 = Polygon(poly2)

get_intersection(polygon1, polygon2)
``````

I get:

``````area per unit??  1.8490791593183024e-06
``````

Is there any quick conversion of whatever I received in square metes?

EDIT The following function takes too long to transform the coordinates.

``````def transform_to_meters(coordinates):
coord_tuple = []
for (x,y) in coordinates:
inProj = Proj('epsg:4326')
outProj = Proj('epsg:3857')
lat, lon = transform(inProj,outProj,x,y)
coord_tuple.append((lat, lon))
return coord_tuple
``````
• You absolutely do not want to use 3857 (Web Mercator) for distance or area. It is always wrong (the poles are infinitely far from the Equator). Many tools can calculate geodesic area from 4326 coordinates; use one of them. Commented Apr 29, 2020 at 11:41
• I cannot use 4326 for the intersection as it provides me degrees, what shall I use then? Commented Apr 29, 2020 at 11:42
• I'd use a tool that does provide geodesic area (PostGIS or ArcPy), or reproject each geometry to a local projection centered on the shape or use a single best-fit CRS for the study area. Commented Apr 29, 2020 at 11:50
• For calculating area you need to use equal area projection since then you will have lowest distortion of area. If you have some local equal area projection that'd be great, otherwise you can use EPSG:3035 or global Albers equal area projection. Commented Apr 29, 2020 at 12:40
• Side note: At least one of your data sets is not in EPSG 4326.
– Erik
Commented Apr 30, 2020 at 6:41

No. One degree from east to west at the equator is about 110 km, while at the polar circle it is about 38 km. Geographical coordinates are not meant for area and distance calculations.

• What would be the fastest way to calculate the area in square meters? The transform in the edit is too slow for me. Commented Apr 29, 2020 at 11:17
• Get the data in any m-based CRS. There is no other way.
– Erik
Commented Apr 29, 2020 at 11:18
• I cannot, this is all I got. Commented Apr 29, 2020 at 11:18
– Erik
Commented Apr 29, 2020 at 11:19
• Personally I'd use QGIS.
– Erik
Commented Apr 29, 2020 at 11:22

If you want to get the area in square meters you need to use a projection which uses meter as a unit. You cannot use WGS84. I propose to use any equal area projection.

I'm not sure why to you transform each point separately and then build coordinates list.

I responded to another post from you where you could take pieces of the code.