Find the area between two curves plotted in matplotlib (fill_between area)

Question

I have a list of x and y values for two curves, both having weird shapes, and I don't have a function for any of them. I need to do two things:

1. Plot it and shade the area between the curves like the image below.
2. Find the total area of this shaded region between the curves.

I'm able to plot and shade the area between those curves with fill_between and fill_betweenx in matplotlib, but I have no idea on how to calculate the exact area between them, specially because I don't have a function for any of those curves.
Any ideas?

I looked everywhere and can't find a simple solution for this. I'm quite desperate, so any help is much appreciated.

Thank you very much!

---

EDIT: For future reference (in case anyone runs into the same problem), here is how I've solved this: connected the first and last node/point of each curve together, resulting in a big weird-shaped polygon, then used shapely to calculate the polygon's area automatically, which is the exact area between the curves, no matter which way they go or how nonlinear they are. Works like a charm! :)

Here is my code:

from shapely.geometry import Polygon

x_y_curve1 = [(0.121,0.232),(2.898,4.554),(7.865,9.987)] #these are your points for curve 1 (I just put some random numbers)
x_y_curve2 = [(1.221,1.232),(3.898,5.554),(8.865,7.987)] #these are your points for curve 2 (I just put some random numbers)

polygon_points = [] #creates a empty list where we will append the points to create the polygon

for xyvalue in x_y_curve1:
    polygon_points.append([xyvalue[0],xyvalue[1]]) #append all xy points for curve 1

for xyvalue in x_y_curve2[::-1]:
    polygon_points.append([xyvalue[0],xyvalue[1]]) #append all xy points for curve 2 in the reverse order (from last point to first point)

for xyvalue in x_y_curve1[0:1]:
    polygon_points.append([xyvalue[0],xyvalue[1]]) #append the first point in curve 1 again, to it "closes" the polygon

polygon = Polygon(polygon_points)
area = polygon.area
print(area)

EDIT 2: Thank you for the answers. Like Kyle explained, this only works for positive values. If your curves go below 0 (which is not my case, as showed in the example chart), then you would have to work with absolute numbers.

Answer 1

The area calculation is straightforward in blocks where the two curves don't intersect: thats the trapezium as has been pointed out above. If they intersect, then you create two triangles between x[i] and x[i+1], and you should add the area of the two. If you want to do it directly, you should handle the two cases separately. Here's a basic working example to solve your problem. First, I will start with some fake data:

#!/usr/bin/python
import numpy as np

# let us generate fake test data
x = np.arange(10)
y1 = np.random.rand(10) * 20
y2 = np.random.rand(10) * 20

Now, the main code. Based on your plot, looks like you have y1 and y2 defined at the same X points. Then we define,

z = y1-y2
dx = x[1:] - x[:-1]
cross_test = np.sign(z[:-1] * z[1:])

cross_test will be negative whenever the two graphs cross. At these points, we want to calculate the x coordinate of the crossover. For simplicity, I will calculate x coordinates of the intersection of all segments of y. For places where the two curves don't intersect, they will be useless values, and we won't use them anywhere. This just keeps the code easier to understand.

Suppose you have z1 and z2 at x1 and x2, then we are solving for x0 such that z = 0:

# (z2 - z1)/(x2 - x1) = (z0 - z1) / (x0 - x1) = -z1/(x0 - x1)
# x0 = x1 - (x2 - x1) / (z2 - z1) * z1
x_intersect = x[:-1] - dx / (z[1:] - z[:-1]) * z[:-1]
dx_intersect = - dx / (z[1:] - z[:-1]) * z[:-1]

Where the curves don't intersect, area is simply given by:

areas_pos = abs(z[:-1] + z[1:]) * 0.5 * dx # signs of both z are same

Where they intersect, we add areas of both triangles:

areas_neg = 0.5 * dx_intersect * abs(z[:-1]) + 0.5 * (dx - dx_intersect) * abs(z[1:])

Now, the area in each block x[i] to x[i+1] is to be selected, for which I use np.where:

areas = np.where(cross_test < 0, areas_neg, areas_pos)
total_area = np.sum(areas)

That is your desired answer. As has been pointed out above, this will get more complicated if the both the y graphs were defined at different x points. If you want to test this, you can simply plot it (in my test case, y range will be -20 to 20)

negatives = np.where(cross_test < 0)
positives = np.where(cross_test >= 0)
plot(x, y1)
plot(x, y2)
plot(x, z)
plt.vlines(x_intersect[negatives], -20, 20)

Answer 2

Define your two curves as functions f and g that are linear by segment, e.g. between x1 and x2, f(x) = f(x1) + ((x-x1)/(x2-x1))*(f(x2)-f(x1)). Define h(x)=abs(g(x)-f(x)). Then use scipy.integrate.quad to integrate h.

That way you don't need to bother about the intersections. It will do the "trapeze summing" suggested by ch41rmn automatically.