这可能是一个非常简单的问题,但似乎我无法看到它.
我有一个顺时针排序的点列表,想要根据
this使用以下函数计算这些点的质心(凸多边形):
和
def calculateCentroid(raLinks,raNodes, links, nodes):
orderedPointsOfLinks = orderClockwise(raLinks,raNodes, links, nodes)
arg1 = 0
arg2 = 0
Xc = 0
Yc = 0
i = 0
for point in orderedPointsOfLinks:
arg1 += point.Y*(orderedPointsOfLinks[i+1 if i+1<len(orderedPointsOfLinks) else 0].X)
arg2 += (orderedPointsOfLinks[i+1 if i+1<len(orderedPointsOfLinks) else 0].Y)*point.X
Xc += (point.X+(orderedPointsOfLinks[i+1 if i+1<len(orderedPointsOfLinks) else 0].X))*(((orderedPointsOfLinks[i+1 if i+1<len(orderedPointsOfLinks) else 0].Y)*point.X)-(point.Y*(orderedPointsOfLinks[i+1 if i+1<len(orderedPointsOfLinks) else 0].X)))
Yc += (point.Y+(orderedPointsOfLinks[i+1 if i+1<len(orderedPointsOfLinks) else 0].Y))*(((orderedPointsOfLinks[i+1 if i+1<len(orderedPointsOfLinks) else 0].Y)*point.X)-(point.Y*(orderedPointsOfLinks[i+1 if i+1<len(orderedPointsOfLinks) else 0].X)))
i+=1
area = (arg1-arg2)*0.5
print area
X = -Xc/(6*area)
Y = -Yc/(6*area)
print X , " ", Y
使用Arcpy计算面积和centorid表示上述函数计算的面积是正确的,但质心是错误的.
Xc和Yc有什么问题我无法修复它?
如果我按以下方式更改for循环,它可以工作:
for point in orderedPointsOfLinks:
y0 = point.Y
x0 = point.X
x1 = orderedPointsOfLinks[i+1 if i+1<len(orderedPointsOfLinks) else 0].X
y1 = orderedPointsOfLinks[i+1 if i+1<len(orderedPointsOfLinks) else 0].Y
a = x0*y1 - x1*y0
area += a
Xc += (x0+x1)*a
Yc += (y0+y1)*a
i+=1
area *= 0.5
print area
X = Xc/(6*area)
Y = Yc/(6*area)
print X , " ", Y
这是检查代码的节点列表:
[(371623.876, 6159668.714),(371625.994, 6159661.094), (371624.319, 6159654.634), (371619.654, 6159649.86), (371614.194, 6159647.819), (371608.401, 6159648.449), (371601.544, 6159652.652), (371598.77, 6159658.058), (371599.318, 6159665.421), (371603.025, 6159671.805), (371611.372, 6159674.882 ), (371619.417, 6159673.065)]
最佳答案
source
尝试:
import numpy
tp = [(371623.876, 6159668.714),(371625.994, 6159661.094), (371624.319, 6159654.634), (371619.654, 6159649.86),\
(371614.194, 6159647.819), (371608.401, 6159648.449), (371601.544, 6159652.652), (371598.77, 6159658.058), \
(371599.318, 6159665.421), (371603.025, 6159671.805), (371611.372, 6159674.882 ), (371619.417, 6159673.065),(371623.876, 6159668.714)]
# cx = sigma (x[i]+x[i+1])*((x[i]*y[i+1]) - (x[i+1]*y[i] ))
# cy = sigma (y[i]+y[i+1])*((x[i]*y[i+1]) - (x[i+1]*y[i] ))
cx = 0
cy = 0
p = numpy.array(tp)
x = p[:, 0]
y = p[:, 1]
a = x[:-1] * y[1:]
b = y[:-1] * x[1:]
cx = x[:-1] + x[1:]
cy = y[:-1] + y[1:]
tp = tp[:-1] #dont need repeat
def area():
tox=0
toy=0
for i in range(len(tp)):
if i+1 == len(tp):
tox += tp[-1][0]*tp[0][1]
else:
tox += tp[i][0]*tp[i+1][1]
for i in range(len(tp)):
if i+1 == len(tp):
toy += tp[-1][1]*tp[0][0]
else:
toy += tp[i][1]*tp[i+1][0]
return abs(tox-toy)*0.5
ar = area()
Cx = abs(numpy.sum(cx * (a - b)) / (6. * ar))
Cy = abs(numpy.sum(cy * (a - b)) / (6. * ar))
print Cx,Cy
警告 !
tp[0] == tp[-1]
所以:第一个和最后一个坐标是相同的值……