軌跡數(shù)據(jù)處理

由于GPS精度以及系統(tǒng)誤差等原因，造成gps軌跡數(shù)據(jù)像狗啃一樣，不是那么規(guī)則，且大多數(shù)點(diǎn)無法落在道路上，因此這篇文章主要是對(duì)GPS軌跡數(shù)據(jù)進(jìn)行處理。

  1. 原始數(shù)據(jù)為csv格式數(shù)據(jù),具體怎么將csv數(shù)據(jù)轉(zhuǎn)化為空間數(shù)據(jù)就不多贅述。以下為gps數(shù)據(jù)表格，此處已簡(jiǎn)化數(shù)據(jù)規(guī)模，表中只有一條軌跡，便于計(jì)算。
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");

以下是數(shù)據(jù)庫表中數(shù)據(jù)。


image

gps軌跡數(shù)據(jù)可視化顯示如下圖：可以看到數(shù)據(jù)很多鋸齒，且未在落在道路還是上。


image

  1. 數(shù)據(jù)處理流程
    (http://upload-images.jianshu.io/upload_images/8667322-f963c6add8a66bc6.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)

  2. 數(shù)據(jù)去重方法，將經(jīng)緯度相同的點(diǎn)去除。此處采用了窗口函數(shù)，簡(jiǎn)化了算法。因?yàn)榇颂幰庖艘恍┢渌姆椒?，一次有一些多余的?shù)據(jù)與參數(shù)，但是目前未實(shí)現(xiàn)，但也不作刪除，說不定，哪天靈感迸發(fā)，想到了解決方案。以下是數(shù)據(jù)去除的腳本。

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,  --窗口函數(shù)上一條記錄
            time, 
            LEAD (TIME) OVER (PARTITION BY DATE ORDER BY TIME) :: TIME AS leadTime,--窗口函數(shù)下一條記錄
            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個(gè)重復(fù)點(diǎn)，為之后的計(jì)算做鋪墊，效果圖如下


image

4.數(shù)據(jù)平滑采用高斯濾波進(jìn)行平滑，直接上代碼：

--平滑軌跡
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點(diǎn)位前后三十秒的數(shù)據(jù)進(jìn)行加權(quán)平滑
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));
--平滑后的數(shù)據(jù)入庫
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();

以下是平滑后的數(shù)據(jù)，可以看到比原始數(shù)據(jù)明顯少了很多鋸齒，效果還是可以的。


image
  1. 最后將平滑后的數(shù)據(jù)，匹配到道路上去，目前之做到簡(jiǎn)單匹配，有些地方匹配結(jié)果是錯(cuò)誤的，但是目前沒有想到好的方法。以下是匹配代碼：
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
    --軌跡點(diǎn)連成線
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米范圍生成緩沖區(qū)，并計(jì)算與緩沖區(qū)有相交的道路
    --點(diǎn)落在道路上
FOR gpsPoint IN 
SELECT * FROM gps_data_smooth loop 
raise notice '正在處理gid:%',   '----------------------------------' ;
 --計(jì)算離Gps點(diǎn)位最近的道路，并計(jì)算垂點(diǎn)，將垂點(diǎn)入庫
 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);
--數(shù)據(jù)入庫
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 ; $$;

以下數(shù)匹配的結(jié)果圖，可以看到能匹配到道路上去，但是由于精度不是那么可靠，切在轉(zhuǎn)彎處的數(shù)據(jù)匹配也是明顯的錯(cuò)誤，但是目前沒找到好的解決方案。（Ps:目前欲采用軌跡的角度和道路的角度進(jìn)行計(jì)算，但是目前效果不佳，持續(xù)改進(jìn)中）


image

總結(jié)：

  • 用PG進(jìn)行數(shù)據(jù)數(shù)據(jù)處理很方便，尤其是在批量處理上，能夠靈活組裝函數(shù)。
  • 很多函數(shù)用的不是很熟，處理函數(shù)的報(bào)錯(cuò)過程中花費(fèi)太多時(shí)間。
  • 算法還較為簡(jiǎn)單，需要持續(xù)改進(jìn)。
  • PG so funny!
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容