Вычислить площадь заданных координат (x, y)
У меня есть набор точек и хотелось бы знать, есть ли функция (ради удобства и, вероятно, скорости), которая может вычислять область, заключенную в набор точек.
например:
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... и они будут многоугольниками без какой-либо определенной функции.
EDIT:
точки, скорее всего, не будут в каком-либо заданном порядке (по часовой стрелке или против часовой стрелки) и могут быть довольно сложными, поскольку они представляют собой набор координат utm из шейп файла под набором границ
Ответы
Ответ 1
Внедрение формулы Shoelace можно осуществить в 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
цикла делает эту функцию ~ 50X быстрее, чем 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.
Ответ 2
Вы можете использовать формулу шнурки > , например
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)]
Это работает только для простых полигонов
Ответ 3
np.roll()
ответ Махди, я пришел к выводу, что большая часть времени была потрачена на выполнение 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()
Ответ 4
В приведенном выше коде есть ошибка, поскольку на каждой итерации она не принимает абсолютных значений. Вышеприведенный код всегда будет возвращать ноль. (Математически, это разница между принятием подписанного участка или клиновым продуктом и фактической площадью
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
Ответ 5
Это намного проще для правильных многоугольников:
import math
def area_polygon(n, s):
return 0.25 * n * s**2 / math.tan(math.pi/n)
так как формула равна ¼ n s2/tan (π/n).
Учитывая количество сторон, n и длину каждой стороны, s
Ответ 6
Основываясь на
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)
Фокус в том, что первая координата также должна быть последней.
Ответ 7
Ответ maxb дает хорошую производительность, но может легко привести к потере точности, когда значения координат или количество точек велики. Это может быть смягчено простым сдвигом координат:
def polygon_area(x,y):
# coordinate shift
x_ = x - x.mean()
y_ = y - y.mean()
# everything else is the same as maxb 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()
из каждого продукта - Он производит сочетание положительных и отрицательных значений в каждом точечном произведении, которое в значительной степени отменяется.
Смещение координат не изменяет общую площадь, оно просто делает расчет более устойчивым в численном отношении.