5.轨迹抽象，关键点序列生成

研究了两篇模型，接下来继续进行轨迹抽象，试一试模型效果。再然后统计预测效果或重新研究预测模型

现在发现了几个问题：

1.转移概率是从状态序列统计出来的，那么回头验证的时候准确率肯定和转移概率一样。就是转移概率低的点肯定预测不到。需要结合其他信息预测

2.用三维矩阵统计出的马尔可夫链统计出的数量比预计的少，很多都是1，没有代表性。二维矩阵又有丢失信息，二者结合应该是个思路，那干脆更多结合好了。

找到一篇论文说低阶马尔可夫链准确率低，高阶马尔可夫链稀疏率高。

现在看起来关键点的选取没什么可研究的，能代表轨迹就行。更有研究价值的在于数据准备和预测模型，看哪边更容易突破吧。我去查预测模型相关的文献去了。

另，附上关键点序列生成的代码和结果。

代码是在之前的基础上添加的。最后会输出一个OID的序列，做成表（这步没有API，所以后面就手动了），连接上极值面要素，显示XY数据，导出成图层，XY转线，

# -*- coding: UTF-8 -*-# Name: CreateFishnet.py# Description: Creates rectangular cells# import system moduleimport arcpyimport timefrom arcpy import envdef GetDistance(lng1, lat1, lng2, lat2):    def rad(num):        return num/180*math.pi    EARTH_RADIUS = 6378137    radLat1 = rad(lat1);    radLat2 = rad(lat2);    a = radLat1 - radLat2;    b = rad(lng1) - rad(lng2);    s = 2 * math.asin(math.sqrt(math.pow(math.sin(a/2),2) +    math.cos(radLat1)*math.cos(radLat2)*math.pow(math.sin(b/2),2)));    s = s * EARTH_RADIUS;    #s = math.round(s * 10000) / 10000;    return s;def GetMinMax(input):    fc=arcpy.SearchCursor(input)    minx=999;    maxx=0;    miny=999;    maxy=0;    for row in fc:        minx=min(minx,row.getValue("LONGITUDE_"));        maxx=max(maxx,row.getValue("LONGITUDE_"));        miny=min(miny,row.getValue("LATITUDE_D"));        maxy=max(maxy,row.getValue("LATITUDE_D"));    minx=minx-minx%0.01;    maxx=maxx-maxx%0.01+0.02;    miny=miny-miny%0.01;    maxy=maxy-maxy%0.01+0.02;    list=[0,0,0,0];    list[0]=minx;    list[1]=maxx;    list[2]=miny;    list[3]=maxy;    return list;    def nearTop(fea,search_radius):    def get(row0):        feature1.reset();        for row1 in feature1:            if(row1[0]==row0):                return row1[1];    def update(row0,str):        updateFea.reset();        for row1 in updateFea:            if(row1[0]==row0):                if(row1[1]==0):                    return;                else:                    row1[1]=str;                    updateFea.updateRow(row1)                    #arcpy.AddMessage("update:");                return;    def test():        n=0;        for row in neartable:            arcpy.AddMessage("process:"+str(row[0]));            n+=1;            if(get(row[0])>=get(row[1])):                update(row[0],1);                update(row[1],0);            if(n==-1):                return;    temp="temp"+str(int(time.time()))    arcpy.AddField_management(fea, "test", "SHORT")    arcpy.GenerateNearTable_analysis(fea, fea, temp, search_radius, 'LOCATION', 'NO_ANGLE', 'ALL','#')    neartable=arcpy.da.SearchCursor(temp,["IN_FID","NEAR_FID"])    feature1=arcpy.da.SearchCursor(fea,["OID@","Join_Count"])    updateFea=arcpy.da.UpdateCursor(fea,["OID@","test"])                    test();    #arcpy.DeleteFeatures_management(temp)    def clear(fea):    updateFea=arcpy.da.UpdateCursor(fea,["OID@","test"])    for row in updateFea:        if row[1]==0:            updateFea.deleteRow()    updateFea.reset()    arcpy.CopyFeatures_management(fea,fea+"_copy")    return fea+"_copy"    def getToken(fea):    searchFea=arcpy.da.SearchCursor(fea,["NEAR_FID"])    token=""    last=0;    for row in searchFea:        if (row[0]!=last):            last=row[0]            token+=str(last)+","    return token# set workspace environmentenv.workspace = "C:\Users\panda\Desktop\轨迹预测\轨迹.gdb"# Set coordinate system of the output fishnetenv.outputCoordinateSystem = arcpy.SpatialReference("WGS 1984")input = arcpy.GetParameterAsText(0)output = arcpy.GetParameterAsText(1)isLabel=arcpy.GetParameterAsText(2)search_radius_text=arcpy.GetParameterAsText(3)label='NO_LABELS'if isLabel:    label='LABELS'#if(!(search_radius_text==''&&search_radius_text=='#'))search_radius=float(search_radius_text)list=GetMinMax(input);temp="temp"+str(int(time.time()))arcpy.CreateFishnet_management(temp, str(list[0])+' '+str(list[2]), str(list[0])+' '+str(list[2]+1), '0.01', '0.01', '0', '0', str(list[1])+' '+str(list[3]), label, '#', 'POLYGON')#arcpy.DeleteFeatures_management(temp)arcpy.SpatialJoin_analysis(temp, input, output, "JOIN_ONE_TO_ONE", "KEEP_COMMON", "#", "INTERSECT")nearTop(output,search_radius)temp="temp"+str(int(time.time()))arcpy.GenerateNearTable_analysis(input, clear(output), temp, search_radius, 'LOCATION', 'NO_ANGLE', 'ALL','#')arcpy.Sort_management(temp,temp+"_Sort",[["IN_FID","ASCENDING"]])token=getToken(temp+"_Sort")arcpy.AddMessage("token:"+token);

这就是结果，标记的是序号，按数字顺序就是轨迹顺序。如果关键点足够多应该抽象出来的轨迹是可以的，刚刚发现一样的代码换一台电脑跑出来的关键点数量居然不一样（所以这张图的关键点比之前实验的少了）。以后解决。现在先集中精力做理论研究。