点数据集描述性空间统计之六——方向分布统计（标准差椭圆）（Standard DeviationalEllipse)原理及python实现

1.原理

在上一篇圆概率误差统计中（详见：【空间分析之五】点数据集圆概率误差统计CEP（Circle Error Probable ) - ganquan78的博客 - CSDN博客），我们实现了对点数据集点X坐标、Y坐标以及距离的标准差计算，求出了以平均中心为圆心，CEP为半径的概率误差圆。这种统计方法简单易行，但是其缺陷在于，无法统计出点数据的方向趋势。在真实的地理环境下，点数据集的方向分布往往跟真实地理实体的走向具有相关性，比如某种流行病的分布方向与附近的河流走向紧密相关，又比如某种类型的犯罪事件分布方向与某条街道走向紧密相关，等等。这些方向特性可以通过对点数据的方向分布统计来计算。
通过测量平差等理论知识我们知道,点数据的距离标准差:
δd=δx2+δy2\delta_d=\sqrt{\delta_x^2+\delta_y^2} δd=δx2+δy2
具有与坐标系统无关的特性，即无论坐标轴如何旋转，距离标准差始终保持不变，无法反应出在每个方位的标准差。
因此，可以通过定义一个标准差椭圆，来描述点数据集分布方向。描述标准差椭圆需要如下几个参数
1.椭圆的中心：为点数据集的平均中心。
2.椭圆的X轴的方向角（以正北为0度，顺时针旋转到X轴的角度）：
3.椭圆的长半轴的长度。坐标旋转后，点数据集在X轴的标准差大小，是所有方向标准差的最大值；
4.椭圆的短半轴的长度。坐标旋转后，点数据集在Y轴的标准差大小，是所有方向标准差的最小值。

2.计算公式

1）计算点数据集的平均中心，根据【空间分析之一】点数据集平均中心统计（Mean Center） - ganquan78的博客 - CSDN博客
可以求出点数据集的平均中心的X坐标和Y坐标：
X‾和Y‾\overline{X}和 \overline{Y}X和Y
2）计算椭圆长轴方向角。
θ=atan(A+BC)\theta=atan (\dfrac{A+B}{C}) θ=atan(CA+B)
A=∑i=1N(xi−X‾)2−∑i=1N(yi−Y‾)2A=\sum_{i=1}^{N}(x_i- \overline{X})^2-\sum_{i=1}^{N}(y_i- \overline{Y})^2A=i=1∑N(xi−X)2−i=1∑N(yi−Y)2
B=A2+4(∑i=1N(xi−X‾)(yi−Y‾))2B=\sqrt{A^2+4(\sum_{i=1}^N(x_i-\overline{X})(y_i-\overline{Y}))^2 }B=A2+4(i=1∑N(xi−X)(yi−Y))2
C=2∑i=1N(xi−X‾)(yi−Y‾)C=2\sum_{i=1}^N(x_i-\overline{X})(y_i-\overline{Y}) C=2i=1∑N(xi−X)(yi−Y)

3）长、短半轴的标准差值，可以用坐标轴旋转的矩阵求解方式计算

δx=2∑i=1N((Xi−X‾)cosθ−(Yi−Y‾)sinθ)2N\delta_x=\sqrt{2} \sqrt{\dfrac{\sum_{i=1}^{N}((X_i-\overline{X})cos\theta-(Y_i-\overline{Y})sin\theta)^2}{N}} δx=2N∑i=1N((Xi−X)cosθ−(Yi−Y)sinθ)2

δx=2∑i=1N((Xi−X‾)sinθ+(Yi−Y‾)conθ)2N\delta_x=\sqrt{2} \sqrt{\dfrac{\sum_{i=1}^{N}((X_i-\overline{X})sin\theta+(Y_i-\overline{Y})con\theta)^2}{N}} δx=2N∑i=1N((Xi−X)sinθ+(Yi−Y)conθ)2

3.源代码

SSDataObject.SSDataObject("c:/General Sample Data/BALTPOP_export.shp")
#  得到图层中要素的总数
mycount = arcpy.management.GetCount("c:/General Sample Data/BALTPOP_export.shp")
# 通过游标读取数据
ftClass = "c:/General Sample Data/BALTPOP_export.shp"
searchCursor = arcpy.da.SearchCursor(ftClass, ["SHAPE@xy"])
# 分别将每个点的X、Y值写入到list列表中，并算出平均中心
listX = []
listY = []
for row in searchCursor:x, y = row[0]listX.append(x)listY.append(y)
# 求出平均中心
MeanCenterX = sum(listX) / len(listX)
MeanCenterY = sum(listY) / len(listY)
# 求出标准误差椭圆的x轴，沿着顺时针，旋转的角度
nNum = len(listX)
i = 0
listXS2 = []
listYS2 = []
listXYS2 = []
while i < nNum:listXS2.append(math.pow((listX[i] - MeanCenterX), 2))listYS2.append(math.pow((listY[i] - MeanCenterY), 2))listXYS2.append((listX[i] - MeanCenterX) * (listY[i] - MeanCenterY))i += 1
SumXS2 = sum(listXS2)
SumYS2 = sum(listYS2)
SumXYS2 = sum(listXYS2)
# 这里的Angletheta是指的y轴沿顺时针方向旋转，与坐标旋转后X轴的夹角（起算方向为Y轴的正北方向）
# 这里Angletheta如果为负数，是指逆时针旋转，如果为正，表示沿着顺时针旋转
Angletheta = math.atan(((SumXS2 - SumYS2) +math.sqrt(math.pow((SumXS2 - SumYS2), 2) + 4 * SumXYS2 * SumXYS2)) /(2 * SumXYS2))
# AnglethetaDegree这里
# 首先将 Angletheta 的值由弧度变为度
# 计算X轴，沿着顺时针，旋转的角度（起算轴为X轴的正东方向）
RotationDegree = 0  # RotationDegegree 为旋转后，X坐标的地理方位角
AnglethetaDegree = Angletheta * 180 / math.pi
if Angletheta < 0:RotationDegree = 180 - math.fabs(AnglethetaDegree)AnglethetaDegree = 90 - math.fabs(AnglethetaDegree)
if Angletheta > 0:AnglethetaDegree = math.fabs(AnglethetaDegree) + 270RotationDegree = AnglethetaDegree# 计算当X轴沿着顺时针旋转AnglethetaDegree后，在新坐标轴TposedX，TposedY方向的方差
# 这里用矩阵表示坐标轴旋转公式，可以很清晰的反应变化的过程
#  cos（theta）-sin（theta）   x
#  sin(theta)   cos(theta)    y
i = 0
listTPosedX2 = []
listTPosedY2 = []
for i in range(nNum):listTPosedX2.append(((listX[i] - MeanCenterX) * math.cos(AnglethetaDegree * math.pi / 180) -(listY[i] - MeanCenterY) * math.sin(AnglethetaDegree * math.pi / 180)) ** 2)listTPosedY2.append(((listX[i] - MeanCenterX) * math.sin(AnglethetaDegree * math.pi / 180) +(listY[i] - MeanCenterY) * math.cos(AnglethetaDegree * math.pi / 180)) ** 2)i += 1#  tposedSx 为在旋转了当X轴沿着顺时针旋转AnglethetaDegree后，在TposedX轴方向的标准差
#  tposedSy 为在旋转了当X轴沿着顺时针旋转AnglethetaDegree后，在TposedY轴方向的标准差
tposedSx = math.sqrt(2 * sum(listTPosedX2) / nNum)
tposedSy = math.sqrt(2 * sum(listTPosedY2) / nNum)# 生成一个面层，该面层有唯一的一个要素就是，该标准差椭圆
out_path = "c:/General Sample Data/"
out_name = "SDEllipse_test.shp"
out_name_ploygon = "SDEllipse_test2.shp"
geometry_type = "POlygon"
template = "study_quads.shp"
has_m = "DISABLED"
has_z = "DISABLED"
env.workspace = out_path
env.overwriteOutput = True
# MClyr = arcpy.CreateFeatureclass_management(out_path, out_name, geometry_type, None, has_m, has_z)
# 生产一个点图层，点围成了一个椭圆
points = []
pa = arcpy.Array()
i = 0
for i in range(720):  # one point per degree, change if you wisht = math.radians(i / 2)  # 度数换算为弧度，t为旋转了坐标轴i度x = tposedSx * math.cos(t)  # 坐标旋转了i度后，长轴的长度y = tposedSy * math.sin(t)  # 坐标旋转了i度后，短轴的长度theta = math.radians(-AnglethetaDegree)rot_x = MeanCenterX + (x * math.cos(theta)) - (y * math.sin(theta))  # rotate/transpose ellipserot_y = MeanCenterY + (x * math.sin(theta) + (y * math.cos(theta)))  # rotate/transpose ellipsepoints.append(arcpy.PointGeometry(arcpy.Point(rot_x, rot_y)))  # save points to listpa.add(arcpy.Point(rot_x, rot_y))
polygons = arcpy.Polygon(pa)
arcpy.CopyFeatures_management(points, "SDEllipse_test.shp")  # 生成一个点图层
arcpy.PointsToLine_management("SDEllipse_test.shp", "SDEllipse_test1.shp")  # 点图层变为面图层
arcpy.CopyFeatures_management(polygons, "SDEllipse_test2.shp")  # 生成一个面图层# 给面图层赋上属性值
# 使用Arcpy增加字段函数，添加定义好的字段
# 在这个表中增加字段
# 首先要定义各个字段，然后用 addfield命令来实现
McXField = ["McenterX", "Double"]  # 该字段为平均中心点X的坐标值
McYField = ["McenterY", "Double"]  # 该字段为平均中心点Y的坐标值
XStdDistField = ["XStdDist", "Double"]  # 该字段为旋转后X坐标的的标准差
YstdDistField = ["YStdDist", "Double"]  # 该字段为旋转后Y坐标的的标准差
RotationField = ["Rotation", "Double"]  # 该字段为X坐标轴，沿着顺时针方向的旋转角度
arcpy.AddField_management(out_name_ploygon, McXField[0], McXField[1])
arcpy.AddField_management(out_name_ploygon, McYField[0], McYField[1])
arcpy.AddField_management(out_name_ploygon, XStdDistField[0], XStdDistField[1])
arcpy.AddField_management(out_name_ploygon, YstdDistField[0], YstdDistField[1])
arcpy.AddField_management(out_name_ploygon, RotationField[0], RotationField[1])
# 得到第一个要素，并赋值
fields = ["McenterX", "McenterY", "XStdDist", "YStdDist", "Rotation"]
with arcpy.da.UpdateCursor(out_name_ploygon, fields) as cursor:for row in cursor:row[0] = MeanCenterXrow[1] = MeanCenterYrow[2] = tposedSxrow[3] = tposedSyrow[4] = RotationDegreecursor.updateRow(row)

【空间统计之六】点数据集方向分布统计（标准差椭圆）相关推荐

白话空间统计之九：方向分布（标准差椭圆）修正版
文章用红色字体标记出来的内容是修正后的内容,感谢四川的杨同学对我以前的错误提出指正. 终于写到我最喜欢的一个的工具(算法)了,方向分布是虾神我接触的第一个空间统计工具,也是每次讲空间统计必须要讲的一个 ...
新版白话空间统计（25）：方向分布（标准差椭圆）
方向分布是虾神最喜欢的一个空间统计工具,也是最简单明了,但是用处很广的一个点模式的分析中,一般会考察如下五种内容: 1.点的疏密,包括点数据的分布探索,是否一致.均匀或者不均匀. 2.点的方位,包括 ...
方向分布（标准差椭圆）
点模式的分析中,一般会考察如下五种内容: 1.点的疏密,包括点数据的分布探索,是否一致.均匀或者不均匀. 2.点的方位,包括点的分布和方向. 3.点的数量:多少(极值和均值). 4.点的大小:代表的含 ...
从特征值特征向量到方向分析（标准差椭圆）
从特征值特征向量到方向分析(标准差椭圆) 一．特征值与特征向量的意义 Ax=λx 几何直观解释为x向量在矩阵A作用下使得x向量方向不变,且拉伸了λ倍. ...
空间统计（三）聚类分布制图
这组工具中包含众所周知的热点分析工具,通过这个工具我们能捕获到大量数据中的热点和冷点,对我们分析问题有很大的帮助.例如,在犯罪分析中,我们可以研究哪些位置犯罪频繁并且聚集,对增设警力有重要的辅助作用. ...
Seaborn系列(三)：分布统计绘图(distribution)
Seaborn系列目录文章目录 1. 分布统计绘图API概述 2. displot单变量分布图(直方图.核密度.累积分布) 2.1 displot函数绘制单变量分布图 2.2 displot直方图k ...
考研：研究生考试(十五天学完)之【数学考试】—《高等数学-上册/下册》、《线性代数与空间解析几何》、《概率与统计》的研究生学霸重点知识点总结之考试内容各科占比及其知识结构重点
考研:研究生考试(十五天学完)之[数学考试]-<高等数学-上册/下册>.<线性代数与空间解析几何>.<概率与统计>的研究生学霸重点知识点总结之考试内容各科占比及其知 ...
2019北京高考分数分布一览表(成绩分布统计)
北京市2019年普通高考成绩查询系统已开通,考生查询到本人高考成绩,接下来就该考虑填志愿啦.考生成绩分布统计可是填好志愿必不可少的利器,考生和家长们千万看仔细,因为排名对填好志愿意义重大! 成绩分布统 ...
excel 区间人数柱状图_『excle 图表区间计数』excel如何把学生成绩区间分布统计为柱状图...
想要用excel做个图表,就像下面这种,最好还能有区间中值我只会做图表你那种不会这种是:选中数据区>插入>>柱状图>二维,就出现了图表,然后右键图表空白处>选择数据,在 ...

【空间统计之六】点数据集方向分布统计（标准差椭圆）

点数据集描述性空间统计之六——方向分布统计（标准差椭圆）（Standard DeviationalEllipse)原理及python实现

1.原理

2.计算公式

3.源代码

【空间统计之六】点数据集方向分布统计（标准差椭圆）相关推荐

最新文章

热门文章