轨迹数据处理

由于GPS精度以及系统误差等原因,造成gps轨迹数据像狗啃一样,不是那么规则,且大多数点无法落在道路上,因此这篇文章主要是对GPS轨迹数据进行处理。

  1. 原始数据为csv格式数据,具体怎么将csv数据转化为空间数据就不多赘述。以下为gps数据表格,此处已简化数据规模,表中只有一条轨迹,便于计算。
CREATE TABLE "public"."gps_data" (
"gid" int4 DEFAULT nextval('gps_data_gid_seq'::regclass) NOT NULL,
"date" date,
"time" varchar(254) COLLATE "default",
"latitude" numeric,
"longitude" numeric,
"altitude" numeric,
"speed" numeric,
"course" int4,
"type" int4,
"distance" numeric,
"essential" int4,
"geom" "public"."geometry",
CONSTRAINT "gps_data_pkey" PRIMARY KEY ("gid")
)
WITH (OIDS=FALSE)
;
ALTER TABLE "public"."gps_data" OWNER TO "postgres";
CREATE INDEX "gps_data_geom_idx" ON "public"."gps_data" USING gist ("geom");

以下是数据库表中数据。


image

gps轨迹数据可视化显示如下图:可以看到数据很多锯齿,且未在落在道路还是上。


image
  1. 数据处理流程
    (http://upload-images.jianshu.io/upload_images/8667322-f963c6add8a66bc6.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)

  2. 数据去重方法,将经纬度相同的点去除。此处采用了窗口函数,简化了算法。因为此处意淫了一些其他的方法,一次有一些多余的数据与参数,但是目前未实现,但也不作删除,说不定,哪天灵感迸发,想到了解决方案。以下是数据去除的脚本。

delete from gps_data_clean;
DO LANGUAGE plpgsql $$ 
DECLARE  
rec record ; 
declare speed float;
BEGIN
for rec in select *,j.prelength/pretime  prespeed, j.nextlength/nexttime nextspeed from (SELECT
    k.gid,
    k.lagtime,
    extract(epoch FROM (K . TIME :: TIME - K .lagtime)) preTime,
    k.time,
    extract(epoch FROM (K .leadtime - K . TIME :: TIME)) nextTime,
    k.leadtime,
    k.laggeom,
    st_distance(st_transform(k.laggeom, 3857), st_transform(k.geom, 3857)) preLength,
    k.geom,
    st_distance(st_transform(k.geom, 3857), st_transform(k.leadgeom, 3857)) nextLength,
    k.leadgeom
FROM
    (
        SELECT
            gid,
            LAG (TIME) OVER (PARTITION BY DATE ORDER BY TIME) :: TIME AS lagTime,  --窗口函数上一条记录
            time, 
            LEAD (TIME) OVER (PARTITION BY DATE ORDER BY TIME) :: TIME AS leadTime,--窗口函数下一条记录
            geom ,
            LAG (geom) OVER (PARTITION BY DATE ORDER BY TIME) AS lagGeom,
            LEAD (geom) OVER (PARTITION BY DATE ORDER BY TIME) AS leadGeom
        FROM
            gps_data
    ) K) j loop
speed:=(rec.nextspeed+rec.prespeed)/2;

if  speed!=0 then 
 --raise notice '正在处理geom:%',rec;
INSERT INTO "public"."gps_data_clean" ("time","geom") VALUES (rec.time,rec.geom);
else 
raise notice '正在处理geom:%',speed;
end if;
end loop;
END ; $$;

从效果看去除了100个重复点,为之后的计算做铺垫,效果图如下


image

4.数据平滑采用高斯滤波进行平滑,直接上代码:

--平滑轨迹
CREATE
OR REPLACE FUNCTION GetSmoothGpsPt () RETURNS void AS $$
DECLARE vSmoothSpan integer;
declare rec  record;
declare tempRec record;
declare Wi float;
declare Wx float;
declare Wy float;
declare Wa float;
declare sumWX float;
declare sumWY float;
declare sumWA float;
declare sumW float;
declare Latitude float;
declare Longitude float;
declare TimeGap integer;
declare angle float;
BEGIN
vSmoothSpan := 30 ; 
for rec in select *,ST_Azimuth(LAG (geom) OVER (PARTITION BY DATE ORDER BY TIME),LEAD (geom) OVER (PARTITION BY DATE ORDER BY TIME))/(2 * pi()) * 360 angle from gps_data_clean order by time loop
sumWX:= 0; sumWY:= 0; sumWA:= 0;sumW:=0;
--高斯滤波,已Gps点位前后三十秒的数据进行加权平滑
for tempRec in select  *,ST_Azimuth(LAG (geom) OVER (PARTITION BY DATE ORDER BY TIME),LEAD (geom) OVER (PARTITION BY DATE ORDER BY TIME))/(2 * pi()) * 360 angle from gps_data_clean t where t.time::time BETWEEN rec.time::time- interval '30 S' and  rec.time::time+ '30 S' loop
--raise notice '正在处理Longitude:%',tempRec.angle;
TimeGap:=extract(epoch FROM (rec.TIME :: TIME - tempRec.TIME :: TIME ));
Wi:=exp((-1) * TimeGap * TimeGap / (2 * vSmoothSpan * vSmoothSpan));
Wx:= Wi * st_x(tempRec.geom);
Wy:= Wi * st_y(tempRec.geom);
Wa:=Wi*coalesce(tempRec.angle,0);
sumWX = sumWX+Wx;
sumWY = sumWY+Wy;
sumWA=sumWA+Wa;

sumW=sumW+ Wi;
end loop;
Longitude:= sumWX / sumW;
Latitude:= sumWY / sumW;
angle:=sumWA/sumW;

--raise notice '正在处理angle:%',sumWA;
--raise notice '正在处理Longitude:%,Latitude:%,angle:%',Longitude,Latitude,tempRec.angle;
--raise notice '正在处理geom:%',st_astext(ST_GeomFromText('POINT('||Longitude||' '||Latitude||')',4326));
--平滑后的数据入库
INSERT INTO "public"."gps_data_smooth" ("date","time","geom","angle") VALUES ( rec.date,rec.time,ST_GeomFromText('POINT('||Longitude||' '||Latitude||')',4326),angle);
end loop; 
END ; $$ LANGUAGE plpgsql;

select GetSmoothGpsPt();

以下是平滑后的数据,可以看到比原始数据明显少了很多锯齿,效果还是可以的。


image
  1. 最后将平滑后的数据,匹配到道路上去,目前之做到简单匹配,有些地方匹配结果是错误的,但是目前没有想到好的方法。以下是匹配代码:
delete from gps_data_modify;
--select t.geom from gps_data t order by t.time
DO LANGUAGE plpgsql $$ 
DECLARE 
tempRec record ; 
road record ; 
gpsPoint record ; 
rec record ; 
geom geometry ; 
geomArr geometry []; 
point_array geometry []; 
lastPoint geometry ;
angle float;
BEGIN
    --轨迹点连成线
FOR rec IN 
SELECT * FROM   gps_data_smooth ORDER BY TIME loop 
        geomArr := ARRAY [ rec.geom ]; 
        point_array := array_append(point_array, rec.geom) ;
END loop ; --以50米范围生成缓冲区,并计算与缓冲区有相交的道路
    --点落在道路上
FOR gpsPoint IN 
SELECT * FROM gps_data_smooth loop 
raise notice '正在处理gid:%',   '----------------------------------' ;
 --计算离Gps点位最近的道路,并计算垂点,将垂点入库
 FOR tempRec IN 
SELECT *    FROM
            (
                SELECT
                    *,ST_Length (
                        st_transform (
                            ST_ShortestLine (gpsPoint.geom, T .geom),
                            3857
                        )
                    ) min_length,
                    ST_ShortestLine (gpsPoint.geom, T .geom) line
                FROM
                    tianjin_road T
                WHERE
                    ST_Intersects (
                        st_transform (
                            ST_Buffer (
                                st_transform (
                                    ST_MakeLine (point_array),
                                    3857
                                ),
                                40
                            ),
                            4326
                        ),
                        T .geom
                    )
            ) K
        ORDER BY
            K .min_length
        LIMIT 1 loop 


    INSERT INTO "public"."tianjin_road_copy" ("length", "geom")
        VALUES
            (
                ST_Length (
                    st_transform (tempRec.line, 3857)
                ),
                tempRec.line
            ) ; 
lastPoint = ST_ClosestPoint (tempRec.geom, gpsPoint.geom) ; 
angle=ST_Azimuth(gpsPoint.geom,lastPoint);
raise notice '正在处理angle:%',abs(angle/(2 * pi()) * 360-gpsPoint.angle);
--数据入库
INSERT INTO "public"."gps_data_modify" ("date", "time", "geom")
        VALUES
            (
                gpsPoint. DATE,
                gpsPoint. TIME,
                lastPoint
            ) ; 
--raise notice '正在处理gid:%',tempRec.min_length ;
        END loop ;
        END loop ;
        END ; $$;

以下数匹配的结果图,可以看到能匹配到道路上去,但是由于精度不是那么可靠,切在转弯处的数据匹配也是明显的错误,但是目前没找到好的解决方案。(Ps:目前欲采用轨迹的角度和道路的角度进行计算,但是目前效果不佳,持续改进中)


image

总结:

  • 用PG进行数据数据处理很方便,尤其是在批量处理上,能够灵活组装函数。
  • 很多函数用的不是很熟,处理函数的报错过程中花费太多时间。
  • 算法还较为简单,需要持续改进。
  • PG so funny!

推荐阅读更多精彩内容

  • //我所经历的大数据平台发展史(三):互联网时代 • 上篇http://www.infoq.com/cn/arti...
    葡萄喃喃呓语阅读 50,298评论 4 201
  • Spring Cloud为开发人员提供了快速构建分布式系统中一些常见模式的工具(例如配置管理,服务发现,断路器,智...
    卡卡罗2017阅读 131,794评论 18 138
  • Android 自定义View的各种姿势1 Activity的显示之ViewRootImpl详解 Activity...
    passiontim阅读 167,463评论 24 707
  • 宝宝7个月了,在长牙牙,牙床难受啊,每天抓各种东西往嘴里塞。外面买的饼干我不放心,不知道加了多少膨松剂和防腐剂呢!...
    庄经纬阅读 806评论 0 2
  • 1.冒泡排序 2.选择排序 3.插入排序 4.希尔排序 5.归并排序 6.快速排序
    MrThorn阅读 896评论 0 0
  • 李天华提着旅行包从广州火车站出来的时候,天上正下着毛毛细雨,但他还是明显感觉到空气的热度,即将进入7月的广州,让李...
    肖燚阅读 326评论 0 0
  • 高考的前几天,考场号投影在了大屏幕,她从试卷堆里抬起头来,慌忙地在抽屉里寻找眼镜盒,戴上眼镜,找到自己的名字。 1...
    6e7abdae38c4阅读 99评论 0 0