Find positive area between two curves in matplotlib

area_between_curves.py

import math
import numpy as np
import matplotlib.pyplot as plt

from shapely.geometry import Polygon, LineString
from shapely.ops import unary_union, polygonize

# I found it very challenging to find code online that properly finds the positive area between two curves.
# Here is an example that works: we want the area between a sin function and y = 0, which for one full period of a sin function is 4.
# This script calculates this area to be 3.994... close enough!
# In order to make this work, the script calculates the area of each separate polygon (created by the intersection of the two functions) separately.

x = np.arange(0, 4 * math.pi, 0.01)
y1 = np.sin(x)

fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)

pc = ax1.fill_between(x, 0, y1)

sum_area = 0

for path in pc.get_paths():
    for vertices in path.to_polygons():
        vertices_array = np.array(vertices)

        # Create a Shapely Polygon from the vertices
        polygon = Polygon(vertices_array)

        # original data
        ls = LineString(np.c_[x, y1])
        # closed, non-simple
        lr = LineString(ls.coords[:] + ls.coords[0:1])
        assert(lr.is_simple == False)
        
        # Create MultiLineString
        mls = unary_union(lr)
        
        # Sum over LineStrings in mls to get total area
        area_cal =[]

        for polygon in polygonize(mls):
            area_cal.append(polygon.area)
            area_poly = (np.asarray(area_cal).sum())
            x,y = polygon.exterior.xy
            ax2.plot(x,y)
            
        sum_area += area_poly

print(sum_area)

plt.show()