Рассчитать площадь полигона по заданным (х, у) координатам
У меня есть набор точек, и я хотел бы знать, есть ли функция (для удобства и, возможно, для скорости), которая может вычислить площадь, ограниченную набором точек.
например:
x = np.arange(0,1,0.001)
y = np.sqrt(1-x**2)
points = zip(x,y)
дано points площадь должна быть примерно равна (pi-2)/4, Может быть, для этого есть что-то от scipy, matplotlib, numpy, shapely и т. Д.? Я не буду сталкиваться с отрицательными значениями для координат x или y... и они будут многоугольниками без какой-либо определенной функции.
РЕДАКТИРОВАТЬ:
Точки, скорее всего, не будут находиться в каком-либо указанном порядке (по часовой стрелке или против часовой стрелки) и могут быть довольно сложными, поскольку они представляют собой набор координат utm из шейп-файла под набором границ.
13 ответов
Внедрение формулы шнурка может быть сделано в Numpy, Предполагая эти вершины:
import numpy as np
x = np.arange(0,1,0.001)
y = np.sqrt(1-x**2)
Мы можем переопределить функцию в numpy, чтобы найти область:
def PolyArea(x,y):
return 0.5*np.abs(np.dot(x,np.roll(y,1))-np.dot(y,np.roll(x,1)))
И получаю результаты:
print PolyArea(x,y)
# 0.26353377782163534
Как избежать for цикл делает эту функцию в 50 раз быстрее, чем PolygonArea:
%timeit PolyArea(x,y)
# 10000 loops, best of 3: 42 µs per loop
%timeit PolygonArea(zip(x,y))
# 100 loops, best of 3: 2.09 ms per loop.
Время сделано в блокноте Jupyter.
Вы можете использовать формулу шнурка, например,
def PolygonArea(corners):
n = len(corners) # of corners
area = 0.0
for i in range(n):
j = (i + 1) % n
area += corners[i][0] * corners[j][1]
area -= corners[j][0] * corners[i][1]
area = abs(area) / 2.0
return area
# examples
corners = [(2.0, 1.0), (4.0, 5.0), (7.0, 8.0)]
Это работает только для простых полигонов
Если у вас есть многоугольник с отверстиями: рассчитайте площадь внешнего кольца и вычтите области внутренних колец
Если у вас есть самопересекающиеся кольца: вы должны разложить их на простые сектора
Анализируя ответ Махди, я пришел к выводу, что большая часть времени была потрачена на np.roll(), Убрав необходимость в броске и все еще используя numpy, я сократил время выполнения до 4-5 мкс на цикл по сравнению с 41 мкс Махди (для сравнения, функция Махди в среднем занимала 37 мкс на моей машине).
def polygon_area(x,y):
correction = x[-1] * y[0] - y[-1]* x[0]
main_area = np.dot(x[:-1], y[1:]) - np.dot(y[:-1], x[1:])
return 0.5*np.abs(main_area + correction)
Вычисляя корректирующий член, а затем нарезая массивы, нет необходимости катиться или создавать новый массив.
тесты:
10000 iterations
PolyArea(x,y): 37.075µs per loop
polygon_area(x,y): 4.665µs per loop
Время было сделано с использованием time модуль и time.clock()
Ответ maxb дает хорошую производительность, но может легко привести к потере точности, когда значения координат или количество точек велики. Это может быть смягчено простым сдвигом координат:
def polygon_area(x,y):
# coordinate shift
x_ = x - x.mean()
y_ = y - y.mean()
# everything else is the same as maxb's code
correction = x_[-1] * y_[0] - y_[-1]* x_[0]
main_area = np.dot(x_[:-1], y_[1:]) - np.dot(y_[:-1], x_[1:])
return 0.5*np.abs(main_area + correction)
Например, общей системой географической привязки является UTM, которая может иметь (x,y) координаты (488685.984, 7133035.984), Произведение этих двух значений 3485814708748.448, Вы можете видеть, что этот единственный продукт уже на грани точности (он имеет то же количество десятичных разрядов, что и входные данные). Добавление лишь нескольких из этих продуктов, не говоря уже о тысячах, приведет к потере точности.
Простой способ смягчить это - сдвинуть многоугольник от больших положительных координат к чему-то ближе к (0,0), например, вычитая центроид, как в коде выше. Это помогает двумя способами:
- Это устраняет фактор
x.mean() * y.mean()с каждого продукта - Он производит сочетание положительных и отрицательных значений в каждом точечном произведении, которое в значительной степени отменяется.
Сдвиг координат не изменяет общую площадь, он просто делает расчет более устойчивым в численном отношении.
Если точки x и y являются целыми числами, cv2.contourArea() в OpenCV предоставляет альтернативный метод.
points = np.array([[0,0],[10,0],[10,10],[0,10]])
area = cv2.contourArea(points)
print(area) # 100.0
, где аргумент (точки) представляет собой массив numpy (его dtype должен быть int), представляющий вершины многоугольника: [[x1,y1],[x2,y2], ...]
В приведенном выше коде есть ошибка, так как он не принимает абсолютные значения на каждой итерации. Приведенный выше код всегда будет возвращать ноль. (Математически, это разница между взятием продукта со знаком или клина и реальной областью http://en.wikipedia.org/wiki/Exterior_algebra.) Вот некоторый альтернативный код.
def area(vertices):
n = len(vertices) # of corners
a = 0.0
for i in range(n):
j = (i + 1) % n
a += abs(vertices[i][0] * vertices[j][1]-vertices[j][0] * vertices[i][1])
result = a / 2.0
return result
Быстрее использовать shapely.geometry.Polygon а не рассчитывать самостоятельно.
from shapely.geometry import Polygon
import numpy as np
def PolyArea(x,y):
return 0.5*np.abs(np.dot(x,np.roll(y,1))-np.dot(y,np.roll(x,1)))
coords = np.random.rand(6, 2)
x, y = coords[:, 0], coords[:, 1]
С этими кодами и делайте %timeit：
%timeit PolyArea(x,y)
46.4 µs ± 2.24 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit Polygon(coords).area
20.2 µs ± 414 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
Я сравнил все предлагаемые здесь решения с результатом метода площади Shapely, они имели правильную целую часть, но десятичные числа различались. Только решение @Trenton дало правильный результат.
Теперь, улучшая ответ @Trenton на координаты процесса в виде списка кортежей, я придумал следующее:
import numpy as np
def polygon_area(coords):
# get x and y in vectors
x = [point[0] for point in coords]
y = [point[1] for point in coords]
# shift coordinates
x_ = x - np.mean(x)
y_ = y - np.mean(y)
# calculate area
correction = x_[-1] * y_[0] - y_[-1] * x_[0]
main_area = np.dot(x_[:-1], y_[1:]) - np.dot(y_[:-1], x_[1:])
return 0.5 * np.abs(main_area + correction)
#### Example output
coords = [(385495.19520441635, 6466826.196947694), (385496.1951836388, 6466826.196947694), (385496.1951836388, 6466825.196929455), (385495.19520441635, 6466825.196929455), (385495.19520441635, 6466826.196947694)]
Shapely's area method: 0.9999974610685296
@Trenton's area method: 0.9999974610685296
Немного поздно, но не рассматривали ли вы просто использование sympy?
простой код:
from sympy import Polygon
a = Polygon((0, 0), (2, 0), (2, 2), (0, 2)).area
print(a)
На основе
https://www.mathsisfun.com/geometry/area-irregular-polygons.html
def _area_(coords):
t=0
for count in range(len(coords)-1):
y = coords[count+1][1] + coords[count][1]
x = coords[count+1][0] - coords[count][0]
z = y * x
t += z
return abs(t/2.0)
a=[(5.09,5.8), (1.68,4.9), (1.48,1.38), (4.76,0.1), (7.0,2.83), (5.09,5.8)]
print _area_(a)
Хитрость в том, что первая координата также должна быть последней.
Это намного проще для правильных многоугольников:
import math
def area_polygon(n, s):
return 0.25 * n * s**2 / math.tan(math.pi/n)
так как формула ¼ n s2 / tan(π/n). Учитывая количество сторон, n, и длину каждой стороны, с
def find_int_coordinates(n: int, coords: list[list[int]]) -> float:
rez = 0
x, y = coords[n - 1]
for coord in coords:
rez += (x + coord[0]) * (y - coord[1])
x, y = coord
return abs(rez / 2)
Формула шнурков основана на вычислении внутренних треугольников на основе последовательных точек вокруг многоугольника. Было бы более информативно основывать код на этом объяснении. Площадь многоугольника — это просто длина векторного произведения (определителя) двух соседних векторов, деленная на два и сложенная.
https://youtu.be/0KjG8Pg6LGk?t=213
import numpy as np, numpy.linalg as lin
def area(pts):
ps = np.array([0.5 * lin.det(np.vstack((pts[i], pts[i+1]))) for i in range(len(pts)-1)])
s = np.sum(ps)
p1,p2 = pts[-1],pts[0] # cycle back, last pt with the first
s += 0.5 * lin.det(np.vstack((p1,p2)))
return np.abs(s)
points = np.array([[0,0],[10,0],[10,10],[0,10]])
area(points) # 100