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

Question

我有一个 x 和 y 的值的列表，两个曲线都有奇怪的形状，并且我没有任何功能。我需要做两件事：

1. 绘制它并像下图一样遮盖曲线之间的区域。
2. 找到曲线之间该阴影区域的总面积。

我能够使用 fill_betweeen 和 fill_betweenx matplotlib 绘制和阴影这些曲线之间的区域，但是我不知道如何计算计算它们之间的确切区域，特别是因为我没有任何这些曲线功能。
有什么想法吗？

我到处都看，找不到一个简单的解决方案。我非常绝望，因此，任何帮助都非常感谢。

非常感谢！

---

编辑：为将来的参考（如果任何人遇到相同的问题），这就是我如何解决此问题：将每个曲线的第一个和最后一个节点/点连接在一起，从而导致一个大怪异形状的多边形，然后使用 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：感谢您的答案。就像凯尔（Kyle）所解释的那样，这仅适用于积极价值。如果您的曲线低于0（如示例图表中所示，这不是我的情况），那么您必须使用绝对数字。

Answer 1

该区域计算在两条曲线没有相交的块中很简单：那是上面指出的梯形。如果它们相交，则在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.：

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)

Answer 2

将两条曲线定义为函数 f 和 g ，它们是按段线性的，例如 x1 和 x2 之间代码> f（x）= f（x1） +（（x-x1）/（x2-x1））*（f（x2）-f（x1））。
定义 h（x）= abs（g（x）-f（x））。然后使用 scipy.integrate.quad 集成h。

这样，您就不需要打扰交叉点了。它将自动执行CH41RMN建议的“空中求和”。

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

我遇到了同样的问题。下面的答案是基于问题作者的尝试。但是，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)

因此，解决方案是根据线相交的位置将多边形分成较小的片段。然后使用for循环添加这些：

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 5

PYPI库中的 aint_between_two_curves 函数

Answer 6

一般多边形区域的直接应用（参见 shoelace公式）简单而快速，矢量化计算：

def area(p):
    # for p: 2D vertices of a polygon:
    # area = 1/2 abs(sum(p0 ^ p1 + p1 ^ p2 + ... + pn-1 ^ p0))
    # where ^ is the cross product
    return np.abs(np.cross(p, np.roll(p, 1, axis=0)).sum()) / 2

在两条曲线之间的区域应用。在此示例中，我们甚至没有匹配的X坐标！

np.random.seed(0)
n0 = 10
n1 = 15
xy0 = np.c_[np.linspace(0, 10, n0), np.random.uniform(0, 10, n0)]
xy1 = np.c_[np.linspace(0, 10, n1), np.random.uniform(0, 10, n1)]

p = np.r_[xy0, xy1[::-1]]
>>> area(p)
4.9786...

图：

plt.plot(*xy0.T, 'b-')
plt.plot(*xy1.T, 'r-')
p = np.r_[xy0, xy1[::-1]]
plt.fill(*p.T, alpha=.2)

速度

两种曲线具有100万点的

n = 1_000_000
xy0 = np.c_[np.linspace(0, 10, n), np.random.uniform(0, 10, n)]
xy1 = np.c_[np.linspace(0, 10, n), np.random.uniform(0, 10, n)]

%timeit area(np.r_[xy0, xy1[::-1]])
# 42.9 ms ± 140 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

：多边形区域计算的简单效果

# say:
p = np.array([[0, 3], [1, 0], [3, 3], [1, 3], [1, 2]])

p_closed = np.r_[p, p[:1]]
fig, axes = plt.subplots(ncols=2, figsize=(10, 5), subplot_kw=dict(box_aspect=1), sharex=True)
ax = axes[0]
ax.set_aspect('equal')
ax.plot(*p_closed.T, '.-')
ax.fill(*p_closed.T, alpha=0.6)
center = p.mean(0)
txtkwargs = dict(ha='center', va='center')
ax.text(*center, f'{area(p):.2f}', **txtkwargs)
ax = axes[1]
ax.set_aspect('equal')
for a, b in zip(p_closed, p_closed[1:]):
    ar = 1/2 * np.cross(a, b)
    pos = ar >= 0
    tri = np.c_[(0,0), a, b, (0,0)].T
    # shrink a bit to make individual triangles easier to visually identify
    center = tri.mean(0)
    tri = (tri - center)*0.95 + center
    c = 'b' if pos else 'r'
    ax.plot(*tri.T, 'k')
    ax.fill(*tri.T, c, alpha=0.2, zorder=2 - pos)
    t = ax.text(*center, f'{ar:.1f}', color=c, fontsize=8, **txtkwargs)
    t.set_bbox(dict(facecolor='white', alpha=0.8, edgecolor='none'))

plt.tight_layout()

//i.sstatic.net/hqyg3.png“ alt =”在此处输入图像描述”>

Answer 7

该线程上的任何答案都没有完成我想要的，因此我将一些不同的资源合并在一起。我明确希望获得两个相交功能的正区域，其中许多方法在取消区域时存在问题。在此示例中，我正在计算犯罪曲线和y = 0之间的区域，该区域为4个。

要点： https://gist.github.com/amstrudy/0bddab089c7ee36ab49aec8d2363648d

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()

Answer 8

一种优雅的方法是将函数和 shapely.demetry 软件包。 fill_betthe（）返回a 对象，您可以从中获得每个多边形的路径。好事是您甚至可以分别计算 y2＆gt; y1 和 y2＆lt; y1 的区域。

import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import Polygon

n = 10
y1 = np.random.uniform(0, 10, n)
y2 = np.random.uniform(0, 10, n)

fig, ax = plt.subplots()
poly_collect = ax.fill_between(y1, y2)
paths = poly_collect.get_paths()
area = 0
for path in paths:
    poly = Polygon(path.vertices)
    area += poly.area