由于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"); 以下是数据库表中数据。 gps轨迹数据可视化显示如下图：可以看到数据很多锯齿，且未在落在道路还是上。 2、数据处理流程 3、数据去重方法，将经纬度相同的点去除。此处采用了窗口函数，简化了算法。因为此处意淫了一些其他的方法，一次有一些多余的数据与参数，但是目前未实现，但也不作删除，说不定，哪天灵感迸发，想到了解决方案。以下是数据去除的脚本。 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个重复点，为之后的计算做铺垫，效果图如下 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(); 以下是平滑后的数据，可以看到比原始数据明显少了很多锯齿，效果还是可以的。 5、最后将平滑后的数据，匹配到道路上去，目前之做到简单匹配，有些地方匹配结果是错误的，但是目前没有想到好的方法。以下是匹配代码： 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:目前欲采用轨迹的角度和道路的角度进行计算，但是目前效果不佳，持续改进中）总结： 1.用PG进行数据数据处理很方便，尤其是在批量处理上，能够灵活组装函数。 2.很多函数用的不是很熟，处理函数的报错过程中花费太多时间。 3.算法还较为简单，需要持续改进。

PostGIS处理轨迹数据

vip888

由于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");