【空间统计之六】点数据集⽅向分布统计（标准差椭圆）
点数据集描述性空间统计之六——⽅向分布统计（标准差椭圆）（Standard
DeviationalEllip)原理及python 实现
1.原理
在上⼀篇圆概率误差统计中（详见：），我们实现了对点数据集点X坐标、Y坐标以及距离的标准差计算，求出了以平均中⼼为圆⼼，CEP 为半径的概率误差圆。这种统计⽅法简单易⾏，但是其缺陷在于，⽆法统计出点数据的⽅向趋势。在真实的地理环境下，点数据集的⽅向分布往往跟真实地理实体的⾛向具有相关性，⽐如某种流⾏病的分布⽅向与附近的河流⾛向紧密相关，⼜⽐如某种类型的犯罪事件分布⽅向与某条街道⾛向紧密相关，等等。这些⽅向特性可以通过对点数据的⽅向分布统计来计算。
现代教育的特点通过测量平差等理论知识我们知道，点数据的距离标准差:
具有与坐标系统⽆关的特性，即⽆论坐标轴如何旋转，距离标准差始终保持不变，⽆法反应出在每个⽅位的标准差。
因此，可以通过定义⼀个标准差椭圆，来描述点数据集分布⽅向。描述标准差椭圆需要如下⼏个参数
1.椭圆的中⼼：为点数据集的平均中⼼。
2.椭圆的X轴的⽅向⾓（以正北为0度，顺时针旋转到X轴的⾓度）：
3.椭圆的长半轴的长度。坐标旋转后，点数据集在X轴的标准差⼤⼩，是所有⽅向标准差的最⼤值；
4.椭圆的短半轴的长度。坐标旋转后，点数据集在Y轴的标准差⼤⼩，是所有⽅向标准差的最⼩值。
2.计算公式
1）计算点数据集的平均中⼼，根据
可以求出点数据集的平均中⼼的X坐标和Y坐标：
2）计算椭圆长轴⽅向⾓。
3）长、短半轴的标准差值，可以⽤坐标轴旋转的矩阵求解⽅式计算
δ=d δ+δx 2y 2
和X Y
θ=atan ()
C A +B
A =(x −i =1∑N i )−X 2(y −i =1∑N
i )Y 2B =A +4((x −)(y −))2i =1∑N
i X i Y 2
C =2(x −i =1∑N
i )(y −X i )Y
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"
archCursor = arcpy.da.SearchCursor(ftClass, ["SHAPE@xy"])
# 分别将每个点的X 、Y 值写⼊到list 列表中，并算出平均中⼼
listX = []
listY = []
for row in archCursor:
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) + 270
RotationDegree = AnglethetaDegree
# 计算当X 轴沿着顺时针旋转AnglethetaDegree 后，在新坐标轴TpodX ，TpodY ⽅向的⽅差
# 这⾥⽤矩阵表⽰坐标轴旋转公式，可以很清晰的反应变化的过程
# cos （theta ）-sin （theta ） x
# sin(theta) cos(theta) y
δ=x 2N ((X −)cosθ−(Y −)sinθ)∑i =1N
i X i Y 2δ=x 2N
((X −)sinθ+(Y −)conθ)∑i =1N i
X i Y 2
# sin(theta) cos(theta) y
i = 0
listTPodX2 = []
listTPodY2 = []
for i in range(nNum):
全面整顿打一地名listTPodX2.append(((listX[i] - MeanCenterX) * s(AnglethetaDegree * math.pi / 180) -
(listY[i] - MeanCenterY) * math.sin(AnglethetaDegree * math.pi / 180)) ** 2)
listTPodY2.append(((listX[i] - MeanCenterX) * math.sin(AnglethetaDegree * math.pi / 180) +
(listY[i] - MeanCenterY) * s(AnglethetaDegree * math.pi / 180)) ** 2)
i += 1
# tpodSx 为在旋转了当X轴沿着顺时针旋转AnglethetaDegree后，在TpodX轴⽅向的标准差
# tpodSy 为在旋转了当X轴沿着顺时针旋转AnglethetaDegree后，在TpodY轴⽅向的标准差
tpodSx = math.sqrt(2 * sum(listTPodX2) / nNum)
tpodSy = math.sqrt(2 * sum(listTPodY2) / nNum)
# ⽣成⼀个⾯层，该⾯层有唯⼀的⼀个要素就是，该标准差椭圆
out_path = "c:/General Sample Data/"
out_name = "SDEllip_test.shp"
out_name_ploygon = "SDEllip_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 wish
t = math.radians(i / 2) # 度数换算为弧度，t为旋转了坐标轴i度
x = tpodSx * s(t) # 坐标旋转了i度后，长轴的长度
y = tpodSy * math.sin(t) # 坐标旋转了i度后，短轴的长度
theta = math.radians(-AnglethetaDegree)
rot_x = MeanCenterX + (x * s(theta)) - (y * math.sin(theta)) # rotate/transpo ellip
rot_y = MeanCenterY + (x * math.sin(theta) + (y * s(theta))) # rotate/transpo ellip
points.append(arcpy.PointGeometry(arcpy.Point(rot_x, rot_y))) # save points to list
pa.add(arcpy.Point(rot_x, rot_y))
polygons = arcpy.Polygon(pa)
arcpy.CopyFeatures_management(points, "SDEllip_test.shp") # ⽣成⼀个点图层
arcpy.PointsToLine_management("SDEllip_test.shp", "SDEllip_test1.shp") # 点图层变为⾯图层
arcpy.CopyFeatures_management(polygons, "SDEllip_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] = MeanCenterX
row[1] = MeanCenterY
row[2] = tpodSx孕妇怀孕
row[2] = tpodSx
row[3] = tpodSy
row[4] = RotationDegree cursor.updateRow(row)
