由于GPS精度以及系统误差等原因，造成gps轨迹数据像狗啃一样，不是那么规则，且大多数点无法落在道路上，因此这篇文章主要是对GPS轨迹数据进行处理。
- 原始数据为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
数据去重方法，将经纬度相同的点去除。此处采用了窗口函数，简化了算法。因为此处意淫了一些其他的方法，一次有一些多余的数据与参数，但是目前未实现，但也不作删除，说不定，哪天灵感迸发，想到了解决方案。以下是数据去除的脚本。
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
- 最后将平滑后的数据，匹配到道路上去，目前之做到简单匹配，有些地方匹配结果是错误的，但是目前没有想到好的方法。以下是匹配代码：
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!
