在matplotlib中找到两条曲线之间的区域（填充区域）

Question

我有一组 x 和 y 值，代表两条形状奇怪的曲线，但我没有它们的函数。我需要做两件事：

1. 绘制这些曲线，并像下面的图片那样给曲线之间的区域上色。
2. 计算这两条曲线之间上色区域的总面积。

我可以用 matplotlib 中的 fill_between 和 fill_betweenx 来绘制和上色，但我不知道怎么计算它们之间的确切面积，特别是因为我没有这些曲线的函数。
有没有什么想法？

我到处找都找不到简单的解决方案，我有点绝望，所以任何帮助都非常感谢。

非常感谢你们！

---

编辑：为了将来参考（如果有人遇到同样的问题），这是我解决这个问题的方法：把每条曲线的第一个和最后一个节点/点连接起来，形成一个大奇形怪状的多边形，然后用 shapely 自动计算这个多边形的面积，这就是曲线之间的确切面积，无论它们怎么走或者多么不规则。效果很好！ :)

这是我的代码：

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)

编辑 2：谢谢大家的回答。正如凯尔所解释的，这只适用于正值。如果你的曲线低于0（在我的例子中没有这种情况，如示例图表所示），那么你就需要使用绝对值。

Answer 1

在pypi库similaritymeasures中，有一个叫area_between_two_curves的函数（2018年发布），可能正好能满足你的需求。我在这边做了一个简单的例子，比较了一个函数和一个常数值之间的面积，结果和Excel的计算结果非常接近（误差在2%以内）。我不太确定为什么没有达到100%的一致，可能是我哪里做错了。不过，这个方法还是值得考虑的。

Answer 2

我也遇到过同样的问题。下面的回答是基于提问者的尝试。不过，shapely 这个库不会直接给出紫色多边形的面积。你需要修改代码，把它拆分成几个小的多边形，然后分别计算每个的面积。最后，把这些面积加起来就行了。

两条线之间的面积

看看下面的这两条线：

示例：两条线

如果你运行下面的代码，你会发现面积是零，因为它会把顺时针方向的面积减去逆时针方向的面积：

from shapely.geometry import Polygon

x_y_curve1 = [(1,1),(2,1),(3,3),(4,3)] #these are your points for curve 1 
x_y_curve2 = [(1,3),(2,3),(3,1),(4,1)] #these are your points for curve 2 

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)

所以解决办法就是根据线的交点把多边形拆分成更小的部分。然后用一个循环把这些小部分的面积加起来：

from shapely.geometry import Polygon

x_y_curve1 = [(1,1),(2,1),(3,3),(4,3)] #these are your points for curve 1 
x_y_curve2 = [(1,3),(2,3),(3,1),(4,1)] #these are your points for curve 2 

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

x,y = polygon.exterior.xy
    # original data
ls = LineString(np.c_[x, y])
    # closed, non-simple
lr = LineString(ls.coords[:] + ls.coords[0:1])
lr.is_simple  # False
mls = unary_union(lr)
mls.geom_type  # MultiLineString'

Area_cal =[]

for polygon in polygonize(mls):
    Area_cal.append(polygon.area)
    Area_poly = (np.asarray(Area_cal).sum())
print(Area_poly)

Answer 3

你的数据集很“不错”，因为这两组数据的x坐标是相同的。因此，你可以用一系列梯形来计算面积。

比如，定义这两个函数为f(x)和g(x)，那么在任意两个相邻的x点之间，你会有四个数据点：

(x1, f(x1))-->(x2, f(x2))
(x1, g(x1))-->(x2, g(x2))

然后，梯形的面积可以用下面的公式计算：

A(x1-->x2) = ( f(x1)-g(x1) + f(x2)-g(x2) ) * (x2-x1)/2                         (1)

不过，有个复杂的地方是，公式（1）只适用于简单连通的区域，也就是说，这个区域内不能有交叉：

|\             |\/|
|_|     vs     |/\|

交叉部分的两侧面积需要分别计算。你需要检查你的数据，找出所有的交点，然后把它们的坐标加入到你的坐标列表中。x的顺序必须保持正确。接着，你可以遍历你的简单连通区域列表，计算出梯形的面积总和。

补充：

出于好奇，如果这两组数据的x坐标不同，你可以改为构建三角形。例如：

.____.
|   / \
|  /   \
| /     \
|/       \
._________.

在构建三角形时要避免重叠，因此你需要再次找到交点，并将它们加入到你的有序列表中。三角形每条边的长度可以用勾股定理计算，而三角形的面积可以用海伦公式计算。

Answer 4

首先，把你的两个曲线定义为函数 f 和 g，这两个函数是分段线性的。比如说，在 x1 和 x2 之间，f(x) 的计算方式是这样的：f(x) = f(x1) + ((x-x1)/(x2-x1))*(f(x2)-f(x1))。也就是说，你可以通过已知的两个点来计算这个区间内的值。

接下来，定义一个新的函数 h(x)=abs(g(x)-f(x))，这个函数表示了 g 和 f 之间的差距，也就是它们的距离。

然后，你可以使用 scipy.integrate.quad 来对 h 进行积分。这样一来，你就不需要担心它们的交点了。这个方法会自动进行“梯形求和”，正如 ch41rmn 所建议的那样。

Answer 5

计算面积的过程在两个曲线不相交的地方很简单：就像上面提到的那样，是个梯形。如果它们相交，那么你需要在x[i]和x[i+1]之间创建两个三角形，然后把这两个三角形的面积加起来。如果你想直接计算，就要分别处理这两种情况。下面是一个基本的示例来解决你的问题。首先，我会用一些假数据开始：

#!/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

接下来是主要的代码。根据你的图表，看来你在相同的X点上定义了y1和y2。然后我们定义：

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

cross_test在两个图形交叉时会是负值。在这些点上，我们想要计算交叉的x坐标。为了简单起见，我会计算所有y段的交点x坐标。在两个曲线不相交的地方，这些值就没用了，我们也不会用到它们。这样做只是为了让代码更容易理解。

假设你在x1和x2处有z1和z2，那么我们要解出x0，使得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]

在曲线不相交的地方，面积简单地由以下公式给出：

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

在它们相交的地方，我们需要把两个三角形的面积加起来：

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

现在，从每个区间x[i]到x[i+1]的面积需要被选择，为此我使用np.where：

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

这就是你想要的答案。如上所述，如果两个y图在不同的x点上定义，这个过程会变得更复杂。如果你想测试一下，可以简单地画出来（在我的测试案例中，y的范围是从-20到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)