How to calculte the area between a set of coordinates and a give axis?

Question

let's say I have an array of 20 values, for example: [ 0. , -2., 3., -4., etc.] And if I plot them with seaborn lineplot with y = 0, y = 1 and y = -1, they would look like this:

The question is, if I want to calculate the area of the triangle on the left top (the one above the green line), how would I do that in Python? Thanks in advance

Answer 1

This solution is for calculating area between your curve (given x-values and y-values) and y=1, you can change the code accordingly to your needs.

Create the arrays and plot them

To do this we first have to create our arrays (x, y) and check if the line connecting two given point intersect y=1. First of all we create the arrays:

import random

x = [i for i in range(20)]
y = [random.randint(-5, 10) for i in range(20)]

To plot these points and naming them we can use matplotlib (or whatever you want). We also add our y=1 line:

import matplotlib.pyplot as plt
import string

plt.plot(x, y)          
plt.axhline(y=1)        # Straight line
alphabet = string.ascii_lowercase       # naming points on plot
for i in range(len(x)):
    plt.scatter(x[i], y[i])
    plt.annotate(alphabet[i], (x[i], y[i]))
plt.show()

Find the intersection point

Now we have to find if y=1 intercept our line and the x-value of the point:

def intersect(y1, y2):
    '''
    Find if the line connecting two points intercept y=1
    '''
    if (y1 <= 1 and y2 >= 1) or (y2 <= 1 and y1 >= 1):
        return True


def slope_intercept(x1,y1,x2,y2):
    '''
    Find the slope and intercept for the given coordinates 
    and find x value if y=1
    '''
    m = (y2 - y1) / (x2 - x1)
    q = y1 - m * x1    
    point = (1-q)/m 
    return point


intercept = []          # x values of intercept point
for i in range(len(x)-1):
    if intersect(y[i], y[i+1]):
        print(f"{alphabet[i]}-{alphabet[i+1]} intersect")       # which segment intersect y=1
        point_x = round(slope_intercept(x[i],y[i],x[i+1],y[i+1]), 2)  #round our result to two decimal places
        print(f"Intercepted point is x={point_x}")
        intercept.append(point_x)

Find the area

Now that we have our intersection point, we can calculate area of the given polygon. For function we would have used Integral, but here we have straight lines, so we can use the trapezoidal rule. First of all we calculate area under our curve, than we subtract the area under y=1 for the given interval in order to get the area of the polygon:

from math import ceil, floor

def area_calc(interval):
    '''
    Find area of the given interval 
    using the following formula:
    ((x2 - x1) * (f(x2) + f(x1)) / 2
    '''
    tot = 0
    for i in range(len(interval)-1):      # Calculate area for every trapezium    
        h = round(interval[i+1] - interval[i], 2)
        if i == 0:
            b2 = abs(y[interval[i+1]])    # We need to calculate distances, so we use absolute values
            a = (h * (1 + b2)) / 2
            print(f"[{h=} * (1 + {b2=})] / 2 = {a}")
        elif i == len(interval) - 2:
            b1 = abs(y[interval[i]])
            a = (h * (1 + b1)) / 2
            print(f"[{h=} * ({b1=} + 1)] / 2 = {a}")
        else:
            b1 = abs(y[interval[i]])
            b2 = abs(y[interval[i+1]])
            a = (h * (b1 + b2)) / 2
            print(f"[{h=} * ({b1=} + {b2=})] / 2 = {a}")
        tot += abs(a)                         # Total area will be the sum of all trapezium's areas
    tot -= 1 * (interval[-1] - interval[0])   # Subtract area under y=1 for the given interval
    return round(tot, 2)


for i in range(len(intercept)-1):
    start = intercept[i]
    finish = intercept[i+1]
    interval = [z for z in range(ceil(start), floor(finish)+1)]     # Interval used for calculate area
    interval.insert(0, start)       # Insert intercept point at the beginning
    interval.append(finish)         # Insert last point at the end
    #area = area_calc(interval)
    area = area_calc(interval)
    print(f"{start=} {finish=}: {interval=}     {area=}\n")

Finally we get what we wanted!

Conclusion

It might seems difficult at the beginning, but we have to remember some basic calculus and analytical geometry. Of course you can use whatever array you want and calculate areas between your curve and any other curve now, just change some parts of the code.

I hope my explaination was clear and the code work as expected Have fun!

Final code:

from math import ceil, floor
import matplotlib.pyplot as plt
import string
import random 

def intersect(y1, y2):
    '''
    Find if the line connecting two points intercept y=1
    '''
    if (y1 <= 1 and y2 >= 1) or (y2 <= 1 and y1 >= 1):
        return True


def slope_intercept(x1,y1,x2,y2):
    '''
    Find the slope and intercept for the given coordinates 
    and find x value if y=1
    '''
    m = (y2 - y1) / (x2 - x1)
    q = y1 - m * x1    
    point = (1-q)/m 
    return point


def area_calc(interval):
    '''
    Find area of the given interval 
    using the following formula:
    ((x2 - x1) * (f(x2) + f(x1)) / 2
    '''
    tot = 0
    for i in range(len(interval)-1):
        h = round(interval[i+1] - interval[i], 2)
        if i == 0:
            b2 = abs(y[interval[i+1]])
            a = (h * (1 + b2)) / 2
            print(f"[{h=} * (1 + {b2=})] / 2 = {a}")
        elif i == len(interval) - 2:
            b1 = abs(y[interval[i]])
            a = (h * (1 + b1)) / 2
            print(f"[{h=} * ({b1=} + 1)] / 2 = {a}")
        else:
            b1 = abs(y[interval[i]])
            b2 = abs(y[interval[i+1]])
            a = (h * (b1 + b2)) / 2
            print(f"[{h=} * ({b1=} + {b2=})] / 2 = {a}")
        tot += abs(a)
    tot -= 1 * (interval[-1] - interval[0])
    return round(tot, 2)


x = [i for i in range(20)]
y = [random.randint(-5, 10) for i in range(20)]

plt.plot(x, y)          
plt.axhline(y=1)        # Straight line
alphabet = string.ascii_lowercase       # naming points on plot
intercept = []          # x values of intercept point
for i in range(len(x)):
    plt.scatter(x[i], y[i])
    plt.annotate(alphabet[i], (x[i], y[i]))
for i in range(len(x)-1):
    if intersect(y[i], y[i+1]):
        print(f"{alphabet[i]}-{alphabet[i+1]} intersect")       # which segment intersect y=1
        point_x = round(slope_intercept(x[i],y[i],x[i+1],y[i+1]), 2)
        print(f"Intercepted point is x={point_x}")
        intercept.append(point_x)

print()

for i in range(len(intercept)-1):
    start = intercept[i]
    finish = intercept[i+1]
    interval = [z for z in range(ceil(start), floor(finish)+1)]     # Interval used for calculate area
    interval.insert(0, start)       # Insert intercept point at the beginning
    interval.append(finish)         # Insert last point at the end
    #area = area_calc(interval)
    area = area_calc(interval)
    print(f"{start=} {finish=}: {interval=}     {area=}\n")

plt.show()

Simulation case Simulation results:

a-b intersect
Intercepted point is x=0.83
c-d intersect
Intercepted point is x=2.62
d-e intersect
Intercepted point is x=3.36
e-f intersect
Intercepted point is x=4.75
g-h intersect
Intercepted point is x=6.5
i-j intersect
Intercepted point is x=8.43
j-k intersect
Intercepted point is x=9.31
l-m intersect
Intercepted point is x=11.33
m-n intersect
Intercepted point is x=12.18
n-o intersect
Intercepted point is x=13.6
o-p intersect
Intercepted point is x=14.55
q-r intersect
Intercepted point is x=16.57
r-s intersect
Intercepted point is x=17.5
s-t intersect
Intercepted point is x=18.75

[h=0.17 * (1 + b2=2)] / 2 = 0.255
[h=1 * (b1=2 + b2=9)] / 2 = 5.5
[h=0.62 * (b1=9 + 1)] / 2 = 3.1
start=0.83 finish=2.62: interval=[0.83, 1, 2, 2.62]     area=7.07

[h=0.38 * (1 + b2=4)] / 2 = 0.95
[h=0.36 * (b1=4 + 1)] / 2 = 0.8999999999999999
start=2.62 finish=3.36: interval=[2.62, 3, 3.36]     area=1.11

[h=0.64 * (1 + b2=10)] / 2 = 3.52
[h=0.75 * (b1=10 + 1)] / 2 = 4.125
start=3.36 finish=4.75: interval=[3.36, 4, 4.75]     area=6.25

[h=0.25 * (1 + b2=2)] / 2 = 0.375
[h=1 * (b1=2 + b2=2)] / 2 = 2.0
[h=0.5 * (b1=2 + 1)] / 2 = 0.75
start=4.75 finish=6.5: interval=[4.75, 5, 6, 6.5]     area=1.38

[h=0.5 * (1 + b2=4)] / 2 = 1.25
[h=1 * (b1=4 + b2=4)] / 2 = 4.0
[h=0.43 * (b1=4 + 1)] / 2 = 1.075
start=6.5 finish=8.43: interval=[6.5, 7, 8, 8.43]     area=4.4

[h=0.57 * (1 + b2=3)] / 2 = 1.14
[h=0.31 * (b1=3 + 1)] / 2 = 0.62
start=8.43 finish=9.31: interval=[8.43, 9, 9.31]     area=0.88

[h=0.69 * (1 + b2=10)] / 2 = 3.795
[h=1 * (b1=10 + b2=2)] / 2 = 6.0
[h=0.33 * (b1=2 + 1)] / 2 = 0.495
start=9.31 finish=11.33: interval=[9.31, 10, 11, 11.33]     area=8.27

[h=0.67 * (1 + b2=1)] / 2 = 0.67
[h=0.18 * (b1=1 + 1)] / 2 = 0.18
start=11.33 finish=12.18: interval=[11.33, 12, 12.18]     area=0.0

[h=0.82 * (1 + b2=10)] / 2 = 4.51
[h=0.6 * (b1=10 + 1)] / 2 = 3.3
start=12.18 finish=13.6: interval=[12.18, 13, 13.6]     area=6.39

[h=0.4 * (1 + b2=5)] / 2 = 1.2000000000000002
[h=0.55 * (b1=5 + 1)] / 2 = 1.6500000000000001
start=13.6 finish=14.55: interval=[13.6, 14, 14.55]     area=1.9

[h=0.45 * (1 + b2=6)] / 2 = 1.575
[h=1 * (b1=6 + b2=9)] / 2 = 7.5
[h=0.57 * (b1=9 + 1)] / 2 = 2.8499999999999996
start=14.55 finish=16.57: interval=[14.55, 15, 16, 16.57]     area=9.9

[h=0.43 * (1 + b2=5)] / 2 = 1.29
[h=0.5 * (b1=5 + 1)] / 2 = 1.5
start=16.57 finish=17.5: interval=[16.57, 17, 17.5]     area=1.86

[h=0.5 * (1 + b2=7)] / 2 = 2.0
[h=0.75 * (b1=7 + 1)] / 2 = 3.0
start=17.5 finish=18.75: interval=[17.5, 18, 18.75]     area=3.75