Skip to content

Instantly share code, notes, and snippets.

@amstrudy
Created February 15, 2024 18:53
Show Gist options
  • Save amstrudy/0bddab089c7ee36ab49aec8d2363648d to your computer and use it in GitHub Desktop.
Save amstrudy/0bddab089c7ee36ab49aec8d2363648d to your computer and use it in GitHub Desktop.
Find positive area between two curves in matplotlib
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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment