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

Question

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

Answer 1

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.

Answer 2

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.

Link to your post

Basically, Shapely is providing already a possibility to easily build geometries, transform them and calculate areas.