方向分布(标准差椭圆)算法

方向分布(标准差椭圆)算法

计算过程

  1. 计算平均/加权平均中心

  2. 计算加全协方差矩阵

    image-20240312180836180

  3. 计算协方差矩阵的特征值与特征向量

    image-20240312180948920

  4. 计算椭圆长轴和短轴

image-20240312181040985

  1. 计算旋转角度

    image-20240312181112722

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

结果

image-20240312182816282

对比ArcGIS ToolBox计算结果

image-20240312182546901

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

Winemonk

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值