2.数据抽象，提取关键点

gps本身返回的数据肯定是一对经纬度，大概每隔一段时间（假设半小时）返回一次，然后画到地图上就是一个点，或者一堆点。按时间顺序连接起来就是一条线型的轨迹

如果一个人每天都走相同的路线，那么线轨迹应该是相似的，但是点轨迹还是乱七八糟，我首先要有一个方案来描述轨迹，不只是过去的轨迹，还有将来的轨迹。因为预测结果不可能精确到经纬度。这就是轨迹抽象

抽象后的轨迹会损失一些细节，但是穷举的数量变少了（如果直接穷举经纬度太可怕了，废话）。就是从轨迹中选取一些关键点来描述轨迹，只要把关键点按时间顺序连接能和原始轨迹相似就对了，当然还有考虑易表达，能输入到预测模型里。

按照计划，预处理应该是清洗噪声和补全丢失，对于我这个急于看到结果（就算结果不好也可以）的人来说就先跳过，我直接拿原始数据做热点分析，找到热点，把几个热点做关键点就行（后期改用其他方法以提高准确率）。然后遍历轨迹点，每个点属于哪个热点，或者都不属于。就得到了一条用热点表达的轨迹时间序列。模型还没开始研究，估计这个序列可以输入给模型。

下面是细节：

1.arcgis的热点分析函数，

arcpy.HotSpots_stats (Input_Feature_Class, Input_Field, Output_Feature_Class, Conceptualization_of_Spatial_Relationships, Distance_Method, Standardization, {Distance_Band_or_Threshold_Distance}, {Self_Potential_Field}, {Weights_Matrix_File})

是按要素类的某个字段的值来计算热点的，但是我想直接计算点密度，抱歉没有，只能先聚类，就是创建渔网，做空间连接（统计网格中点的数量），根据JoinCount做热点分析。我都有JoinCount了还要热点分析干啥，直接count越大越热呗。

arcpy.CreateFishnet_management(...)
arcpy.SpatialJoin_analysis(...)

但是我不能只用最热的点表达轨迹啊，那样只能表达出常走的路线，走的少就不表达了？当然要继续思考，我使用极值点表达，就是某个点比周围的点count大就行。二维空间上的极值点（哦，用渔网聚类以后就是网格面了，不是点）。没有提供就自己写。

2.标记极值点（面）

先建立邻近表（找到一定范围内的点），遍历临近点，两点比较JoinCount大的标记1，小的标记0，已经标记0的不再更改。最后标记1的就是极大值点，代码如下

# -*- coding: UTF-8 -*-
# Name: CreateFishnet.py
# Description: Creates rectangular cells

# import system module
import arcpy
import time
from arcpy import env


def 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 near(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,["[e-mail]","Join_Count"])
    updateFea=arcpy.da.UpdateCursor(fea,["[e-mail]","test"])                
    test();
    #arcpy.DeleteFeatures_management(temp)
    
# set workspace environment
env.workspace = "C:\Users\panda\Desktop\轨迹预测\轨迹.gdb"

# Set coordinate system of the output fishnet
env.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")
near(output,search_radius)

3.最后的几个点就是用这种方法获取的关键点了，按轨迹时间顺序抽象成关键点序列就行了，代码还没写，以后再说