方向分布(标准差椭圆)算法
计算过程
-
计算平均/加权平均中心
-
计算加全协方差矩阵
-
计算协方差矩阵的特征值与特征向量
-
计算椭圆长轴和短轴
-
计算旋转角度
Python实现
引库
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
from scipy.stats import chi2
二维点测试数据
data = np.array([
[120.096888 ,31.45833657],
[120.0942437 ,31.45833657],
[120.1021768 ,31.45833657],
[120.0995324 ,31.45833657],
[120.1074656 ,31.45833657],
[120.1048212 ,31.45833657],
[120.1127544 ,31.45833657],
[120.11011 ,31.45833657],
[120.1180432 ,31.45833657],
[120.1153988 ,31.45833657],
[120.123332 ,31.45833657],
[120.1206876 ,31.45833657],
[120.1259763 ,31.45833657],
[120.1021768 ,31.46054801],
[120.0995324 ,31.46054801],
[120.1074656 ,31.46054801],
[120.1048212 ,31.46054801],
[120.1127544 ,31.46054801],
[120.11011 ,31.46054801],
[120.1180432 ,31.46054801],
[120.1153988 ,31.46054801],
[120.123332 ,31.46054801],
[120.1206876 ,31.46054801],
[120.1259763 ,31.46054801],
[120.1074656 ,31.46275945],
[120.1048212 ,31.46275945],
[120.1127544 ,31.46275945],
[120.11011 ,31.46275945],
[120.1180432 ,31.46275945],
[120.1153988 ,31.46275945],
[120.123332 ,31.46275945],
[120.1206876 ,31.46275945],
[120.1259763 ,31.46275945],
[120.1074656 ,31.46497089],
[120.1048212 ,31.46497089],
[120.1127544 ,31.46497089],
[120.11011 ,31.46497089],
[120.1180432 ,31.46497089],
[120.1153988 ,31.46497089],
[120.123332 ,31.46497089],
[120.1206876 ,31.46497089],
[120.1286207 ,31.46497089],
[120.1259763 ,31.46497089],
[120.1312651 ,31.46497089],
[120.1074656 ,31.46718234],
[120.1048212 ,31.46718234],
[120.1127544 ,31.46718234],
[120.11011 ,31.46718234],
[120.1180432 ,31.46718234],
[120.1153988 ,31.46718234],
[120.123332 ,31.46718234],
[120.1206876 ,31.46718234],
[120.1286207 ,31.46718234],
[120.1259763 ,31.46718234],
[120.1312651 ,31.46718234],
[120.1762198 ,31.46718234],
[120.1021768 ,31.46939378],
[120.1074656 ,31.46939378],
[120.1048212 ,31.46939378],
[120.1127544 ,31.46939378],
[120.11011 ,31.46939378],
[120.1180432 ,31.46939378],
[120.1153988 ,31.46939378],
[120.123332 ,31.46939378],
[120.1206876 ,31.46939378],
[120.1286207 ,31.46939378],
[120.1312651 ,31.46939378],
[120.1074656 ,31.47160522],
[120.1048212 ,31.47160522],
[120.1127544 ,31.47160522],
[120.11011 ,31.47160522],
[120.1180432 ,31.47160522],
[120.1153988 ,31.47160522],
[120.123332 ,31.47160522],
[120.1206876 ,31.47160522],
[120.1259763 ,31.47160522],
[120.1074656 ,31.47381667],
[120.1048212 ,31.47381667],
[120.1127544 ,31.47381667],
[120.11011 ,31.47381667],
[120.1180432 ,31.47381667],
[120.1153988 ,31.47381667],
[120.123332 ,31.47381667],
[120.1206876 ,31.47381667],
[120.1259763 ,31.47381667],
[120.1312651 ,31.47381667],
[120.1629978 ,31.47381667],
[120.1762198 ,31.47381667],
[120.1021768 ,31.47602811],
[120.1074656 ,31.47602811],
[120.11011 ,31.47602811],
[120.1180432 ,31.47602811],
[120.1153988 ,31.47602811],
[120.123332 ,31.47602811],
[120.1259763 ,31.47602811],
[120.1312651 ,31.47602811],
[120.1762198 ,31.47602811],
[120.1127544 ,31.47823955],
[120.11011 ,31.47823955],
[120.1180432 ,31.47823955],
[120.1153988 ,31.47823955],
[120.123332 ,31.47823955],
[120.1206876 ,31.47823955],
[120.1259763 ,31.47823955],
[120.1471315 ,31.47823955],
[120.1127544 ,31.480451],
[120.11011 ,31.480451],
[120.1180432 ,31.480451],
[120.1153988 ,31.480451],
[120.123332 ,31.480451],
[120.1206876 ,31.480451],
[120.1127544 ,31.48266244],
[120.11011 ,31.48266244],
[120.1180432 ,31.48266244],
[120.1153988 ,31.48266244],
[120.123332 ,31.48266244],
[120.1074656 ,31.48487388],
[120.1127544 ,31.48487388],
[120.11011 ,31.48487388],
[120.1206876 ,31.48487388],
[120.1259763 ,31.48487388],
[120.1074656 ,31.48708532],
[120.1127544 ,31.48708532],
[120.11011 ,31.48708532],
[120.1180432 ,31.48708532],
[120.1153988 ,31.48708532],
[120.1127544 ,31.48929676],
[120.11011 ,31.48929676],
[120.1180432 ,31.48929676],
[120.1153988 ,31.48929676],
[120.1788641 ,31.48929676],
[120.1894417 ,31.48929676],
[120.1127544 ,31.49150821],
[120.11011 ,31.49150821],
[120.1180432 ,31.49150821],
[120.1153988 ,31.49150821],
[120.1206876 ,31.49150821],
[120.1127544 ,31.49371965],
[120.1180432 ,31.49371965],
[120.1153988 ,31.49371965],
[120.1286207 ,31.49371965],
[120.1180432 ,31.49593109],
[120.1153988 ,31.49593109],
[120.123332 ,31.49593109],
[120.1153988 ,31.49814254],
[120.1206876 ,31.49814254],
[120.157709 ,31.49814254],
[120.1788641 ,31.49814254],
[120.1180432 ,31.50035398],
[120.1206876 ,31.50035398],
[120.1259763 ,31.50035398],
[120.1206876 ,31.50256542],
[120.1418427 ,31.50477687],
[120.123332 ,31.50698831],
[120.1391983 ,31.50698831],
[120.157709 ,31.50698831],
[120.1656422 ,31.50698831],
[120.1629978 ,31.50698831],
[120.1180432 ,31.50919975],
[120.1153988 ,31.50919975],
[120.123332 ,31.50919975],
[120.1206876 ,31.50919975],
[120.1259763 ,31.50919975],
[120.1656422 ,31.50919975],
[120.1180432 ,31.5114112],
[120.1153988 ,31.5114112],
[120.123332 ,31.5114112],
[120.1286207 ,31.5114112],
[120.1259763 ,31.5114112],
[120.1312651 ,31.5114112],
[120.1180432 ,31.51362264],
[120.1153988 ,31.51362264],
[120.1206876 ,31.51362264],
[120.1206876 ,31.51583408],
[120.1286207 ,31.52689129],
[120.1312651 ,31.52689129],
[120.0334227 ,31.32786144],
[120.0545778 ,31.32786144],
[120.00169 ,31.33007288],
[120.0704441 ,31.33007288],
[120.0730885 ,31.36324453],
[120.1127544 ,31.36324453],
[120.0836661 ,31.36766741],
[120.0598666 ,31.38093607],
[120.0651554 ,31.38093607],
[120.0598666 ,31.38314751],
[120.062511 ,31.38314751],
[120.0942437 ,31.38535895],
[120.1391983 ,31.38535895],
[120.0598666 ,31.3875704],
[120.0757329 ,31.38978184],
[120.0889549 ,31.38978184],
[120.1365539 ,31.39199328],
[120.0730885 ,31.39420473],
[120.0757329 ,31.39420473],
[120.0651554 ,31.39641617],
[120.0677998 ,31.39641617],
[120.0704441 ,31.39641617],
[120.0730885 ,31.39641617],
[120.0757329 ,31.39641617],
[120.062511 ,31.39862761],
[120.0677998 ,31.39862761],
[120.0704441 ,31.39862761],
[120.062511 ,31.40083906],
[120.0704441 ,31.40083906],
[120.0730885 ,31.40083906],
[120.0757329 ,31.40083906],
[120.0704441 ,31.4030505],
[120.1365539 ,31.40526194],
[120.062511 ,31.40747338],
[120.1180432 ,31.40747338],
[120.014912 ,31.40968482],
[120.0704441 ,31.40968482],
[120.0730885 ,31.41189627],
[120.0757329 ,31.41189627],
[120.0836661 ,31.41189627],
[120.0889549 ,31.41189627],
[120.1180432 ,31.41189627],
[120.123332 ,31.41189627],
[120.1259763 ,31.41189627],
[120.0413559 ,31.41410771],
[120.0651554 ,31.41410771],
[120.0730885 ,31.41410771],
[120.0757329 ,31.41410771],
[120.0836661 ,31.41410771],
[120.0889549 ,31.41410771],
[120.0915993 ,31.41410771],
[120.1153988 ,31.41410771],
[120.1180432 ,31.41410771],
[120.1206876 ,31.41410771],
[120.0651554 ,31.41631915],
[120.0677998 ,31.41631915],
[120.0836661 ,31.41631915],
[120.0863105 ,31.41631915],
[120.0889549 ,31.41631915],
[120.0915993 ,31.41631915],
[120.062511 ,31.4185306],
[120.0651554 ,31.4185306],
[120.0677998 ,31.4185306],
[120.0863105 ,31.4185306],
[120.0889549 ,31.4185306],
[120.0202007 ,31.42074204],
[120.0228451 ,31.42074204],
[120.0254895 ,31.42074204],
[120.0281339 ,31.42074204],
[120.0651554 ,31.42074204],
[120.0677998 ,31.42074204],
[120.0704441 ,31.42074204],
[120.0863105 ,31.42074204],
[120.0889549 ,31.42074204],
[120.0995324 ,31.42074204],
[120.1021768 ,31.42074204],
[120.0175563 ,31.42295348],
[120.0202007 ,31.42295348],
[120.0228451 ,31.42295348],
[120.0254895 ,31.42295348],
[120.0281339 ,31.42295348],
[120.0651554 ,31.42295348],
[120.0677998 ,31.42295348],
[120.0704441 ,31.42295348],
[120.0730885 ,31.42295348],
[120.0836661 ,31.42295348],
[120.0863105 ,31.42295348],
[120.0889549 ,31.42295348],
[120.0915993 ,31.42295348],
[120.0942437 ,31.42295348],
[120.096888 ,31.42295348],
[120.0995324 ,31.42295348],
[120.1021768 ,31.42295348],
[120.1286207 ,31.42295348],
[120.0175563 ,31.42516493],
[120.0202007 ,31.42516493],
[120.0228451 ,31.42516493],
[120.0254895 ,31.42516493],
[120.0281339 ,31.42516493],
[120.0651554 ,31.42516493],
[120.0677998 ,31.42516493],
[120.0704441 ,31.42516493],
[120.0836661 ,31.42516493],
[120.0863105 ,31.42516493],
[120.0889549 ,31.42516493],
[120.0915993 ,31.42516493],
[120.0942437 ,31.42516493],
[120.096888 ,31.42516493],
[120.0995324 ,31.42516493],
[120.1259763 ,31.42516493],
[120.1286207 ,31.42516493],
[120.0175563 ,31.42737637],
[120.0202007 ,31.42737637],
[120.0228451 ,31.42737637],
[120.0651554 ,31.42737637],
[120.0704441 ,31.42737637],
[120.0730885 ,31.42737637],
[120.0783773 ,31.42737637],
[120.0810217 ,31.42737637],
[120.0836661 ,31.42737637],
[120.0863105 ,31.42737637],
[120.0889549 ,31.42737637],
[120.0915993 ,31.42737637],
[120.0942437 ,31.42737637],
[120.0175563 ,31.42958781],
[120.0202007 ,31.42958781],
[120.0228451 ,31.42958781],
[120.0783773 ,31.42958781],
[120.0810217 ,31.42958781],
[120.0836661 ,31.42958781],
[120.0863105 ,31.42958781],
[120.0889549 ,31.42958781],
[120.0915993 ,31.42958781],
[120.0942437 ,31.42958781],
[120.0519334 ,31.43179925],
[120.062511 ,31.43179925],
[120.0730885 ,31.43179925],
[120.0783773 ,31.43179925],
[120.0836661 ,31.43179925],
[120.0863105 ,31.43179925],
[120.0889549 ,31.43179925],
[120.0915993 ,31.43179925],
[120.0202007 ,31.4340107],
[120.0863105 ,31.4340107],
[120.0889549 ,31.4340107],
[120.0915993 ,31.4340107],
[120.0942437 ,31.4340107],
[120.096888 ,31.4340107],
[120.1153988 ,31.4340107],
[120.1339095 ,31.4340107],
[120.062511 ,31.43622214],
[120.0889549 ,31.43622214],
[120.0915993 ,31.43622214],
[120.1180432 ,31.43622214],
[120.1206876 ,31.43622214],
[120.1259763 ,31.43622214],
[120.1339095 ,31.43622214],
[120.0757329 ,31.43843358],
[120.0889549 ,31.43843358],
[120.0915993 ,31.43843358],
[120.0942437 ,31.43843358],
[120.1206876 ,31.43843358],
[120.123332 ,31.43843358],
[120.1259763 ,31.43843358],
[120.0915993 ,31.44064502],
[120.0942437 ,31.44064502],
[120.123332 ,31.44064502],
[120.1259763 ,31.44064502],
[120.1286207 ,31.44064502],
[120.1312651 ,31.44064502],
[120.1629978 ,31.44064502],
[120.0281339 ,31.44285647],
[120.0889549 ,31.44285647],
[120.0915993 ,31.44285647],
[120.0942437 ,31.44285647],
[120.096888 ,31.44285647],
[120.0995324 ,31.44285647],
[120.1021768 ,31.44285647],
[120.123332 ,31.44285647],
[120.1259763 ,31.44285647],
[120.1286207 ,31.44285647],
[120.1841529 ,31.44285647],
[120.0836661 ,31.44506791],
[120.0889549 ,31.44506791],
[120.0915993 ,31.44506791],
[120.0942437 ,31.44506791],
[120.096888 ,31.44506791],
[120.0995324 ,31.44506791],
[120.1021768 ,31.44506791],
[120.1048212 ,31.44506791],
[120.1180432 ,31.44506791],
[120.1206876 ,31.44506791],
[120.1259763 ,31.44506791],
[120.1286207 ,31.44506791],
[120.1682866 ,31.44506791],
[120.0810217 ,31.44727935],
[120.0836661 ,31.44727935],
[120.0863105 ,31.44727935],
[120.0889549 ,31.44727935],
[120.0915993 ,31.44727935],
[120.0942437 ,31.44727935],
[120.096888 ,31.44727935],
[120.0995324 ,31.44727935],
[120.1021768 ,31.44727935],
[120.1048212 ,31.44727935],
[120.1418427 ,31.44727935],
[120.0889549 ,31.4494908],
[120.0915993 ,31.4494908],
[120.0942437 ,31.4494908],
[120.096888 ,31.4494908],
[120.0995324 ,31.4494908],
[120.1021768 ,31.4494908],
[120.1048212 ,31.4494908],
[120.1074656 ,31.4494908],
[120.0836661 ,31.45170224],
[120.0863105 ,31.45170224],
[120.0889549 ,31.45170224],
[120.0915993 ,31.45170224],
[120.0942437 ,31.45170224],
[120.096888 ,31.45170224],
[120.0995324 ,31.45170224],
[120.1021768 ,31.45170224],
[120.1048212 ,31.45170224],
[120.1074656 ,31.45170224],
[120.11011 ,31.45170224],
[120.1127544 ,31.45170224],
[120.1153988 ,31.45170224],
[120.123332 ,31.45170224],
[120.0836661 ,31.45391368],
[120.0863105 ,31.45391368],
[120.0889549 ,31.45391368],
[120.0915993 ,31.45391368],
[120.0942437 ,31.45391368],
[120.096888 ,31.45391368],
[120.0995324 ,31.45391368],
[120.1021768 ,31.45391368],
[120.1048212 ,31.45391368],
[120.1074656 ,31.45391368],
[120.11011 ,31.45391368],
[120.1127544 ,31.45391368],
[120.1153988 ,31.45391368],
[120.1259763 ,31.45391368],
[120.0836661 ,31.45612513],
[120.0889549 ,31.45612513],
[120.0915993 ,31.45612513],
[120.0942437 ,31.45612513],
[120.096888 ,31.45612513],
[120.0995324 ,31.45612513],
[120.1021768 ,31.45612513],
[120.1048212 ,31.45612513],
[120.1074656 ,31.45612513],
[120.11011 ,31.45612513],
[120.1127544 ,31.45612513],
[120.1153988 ,31.45612513],
[120.1180432 ,31.45612513],
[120.1206876 ,31.45612513],
[120.123332 ,31.45612513],
[120.1259763 ,31.45612513]])
权重
weights = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 100, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 10, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 10, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 10, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 10, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 10, 1, 1, 1, 10, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 10, 10, 10, 10, 10, 10, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1])
标准差椭圆算法
def calculate_ellipse(da, wg, confidence_level):
# 计算均值
mean = np.average(da, axis=0, weights=wg)
# 计算加权协方差矩阵
cov_matrix = np.cov(da, rowvar=False, aweights=wg)
# 协方差矩阵的特征值和特征向量
eigenvalues, eigenvectors = np.linalg.eig(cov_matrix)
# 计算椭圆长短半轴
alpha = np.sqrt(eigenvalues[0] * chi2.ppf(confidence_level, df=2))
beta = np.sqrt(eigenvalues[1] * chi2.ppf(confidence_level, df=2))
# 计算椭圆角度
angle = np.degrees(np.arctan2(*eigenvectors[:, 0][::-1]))
# 绘制椭圆
ellipse = Ellipse(xy=mean, width=alpha * 2, height=beta * 2, angle=angle, fc='None', lw=2)
return ellipse
绘制各个置信度的标准差椭圆
# 绘制图形
plt.figure(figsize=(12, 12))
plt.scatter(data[:, 0], data[:, 1], label='Points')
ellipse_63 = calculate_ellipse(data, weights, 0.63)
ellipse_63.set_edgecolor('red')
ellipse_63.set_label('63%')
plt.gca().add_patch(ellipse_63)
ellipse_90 = calculate_ellipse(data, weights, 0.90)
ellipse_90.set_edgecolor('orange')
ellipse_90.set_label('90%')
plt.gca().add_patch(ellipse_90)
ellipse_98 = calculate_ellipse(data, weights, 0.98)
ellipse_98.set_edgecolor('yellow')
ellipse_98.set_label('98%')
plt.gca().add_patch(ellipse_98)
ellipse_99 = calculate_ellipse(data, weights, 0.9999)
ellipse_99.set_edgecolor('green')
ellipse_99.set_label('99%')
plt.gca().add_patch(ellipse_99)
# 添加图例
plt.legend()
# 设置图形标签和标题
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Confidence Ellipses')
plt.grid(True)
plt.show()