背景

原始的轨迹数据是一个个单独的GPS点，我们的常见需求一般是

把一系列点组合成一条完整的轨迹；
把整条轨迹按条件（时间，速度，停车点等）分割成多个段；
求一段轨迹的统计信息，如起始点，距离，时间等；

备注：阿里云Ganos提供了更完善的轨迹模型和分析功能，但是仅仅用PostGIS原生功能也足以满足我们的基础需求。

使用PostGIS函数

使用到的PostGIS函数如下

数据来源

T-Drive 公开数据源

地点：中国北京，频率：180秒，持续时长：7天，区域大小：750平方千米，车辆类型：出租车，车辆数目：10357

数据格式如下

1,2008-02-02 15:36:08,116.51172,39.92123
1,2008-02-02 15:46:08,116.51135,39.93883
1,2008-02-02 15:56:08,116.51627,39.91034

字段含义
taxi id, date time, longitude, latitude

原始数据入库

使用Kettle把数据批量入库，表格式如下

drop table if exists demo.t_taxi_gps cascade;
create table demo.t_taxi_gps (tid int, ts timestamp, lng float, lat float);

查询样本数据

select tg.* from demo.t_taxi_gps tg limit 2 ;
  tid  |         ts          |    lng    |   lat
-------+---------------------+-----------+----------
 10012 | 2008-02-02 13:33:43 | 116.46235 | 39.88151
 10012 | 2008-02-02 13:34:12 | 116.46232 | 39.88147
(2 rows)

轨迹分段

业务需求是，根据车辆ID+日期对历史轨迹进行分段。

PostGIS没有提供原生的轨迹(trajectory)数据类型，常规做法是使用LinestringM来进行模拟。思路如下

使用geometry类型，并设置SRID=4326 （WGS84坐标系）。
使用PointM类型保存具体的轨迹点，其中M保存时间戳的unix格式。
使用LinestringM保存一系列点，构造成轨迹数据。

建立轨迹表

drop table if exists demo.t_taxi_trajectory cascade;

-- 使用geometry类型，并设置SRID=4326
create table if not exists demo.t_taxi_trajectory
(
	tid int, 
	dt date, 
	traj geometry(linestringm, 4326),
	primary key (tid, dt)
);

轨迹分段算法

算法很简单，但是有几个需要注意的点

根据车辆+时间戳进行数据去重，目的是保障最后得到合法的轨迹数据。
原始数据中使用了Timestamp Without Time Zone类型保存时间戳，但是转换为Unix时间戳的时候，一定要转换类型带上时区信息，否则会默认为UTC时间。

with
-- 去重并构造 pointM
-- epoch 转换为Unix时间戳要转换类型带上时区信息，否则会默认为UTC时间
point_m as
(
	select distinct on (tid, ts) 
			gt.tid, gt.ts, gt.lng, gt.lat, 
			date_trunc('day', gt.ts)::date dt,
			st_setsrid( st_makepointm( gt.lng, gt.lat, extract(epoch from gt.ts::timestamp with time zone)), 4326) pt
	from demo.t_taxi_gps gt
	order by tid, ts
),
-- 构造轨迹 trajectory
-- array_agg 并带 order by
traj_raw as
(
	select pm.tid, pm.dt,
		st_makeline( array_agg(pm.pt order by pm.ts )::geometry[] ) as traj 	
	from point_m pm 
	group by pm.tid, pm.dt
)
-- 
insert into demo.t_taxi_trajectory
	select tid, dt, traj from traj_raw tr
;

查询轨迹基本信息

查询轨迹的一些基本信息。

注意： 我们把LinestringM保存在 geometry 类型中，如果直接使用 ST_Length() 求轨迹的长度，返回值的单位是度而不是米。所以需要进行类型转换，使用 geography 类型，返回结果单位才是米。

-- 转换为geography类型来计算距离，返回单位是米
select tr.tid, tr.dt,
	st_npoints(tr.traj) as point_count,
	(st_length(tr.traj::geography)/1000)::int as length_kms,
	to_timestamp( st_m( st_startpoint(tr.traj) ) ) ts_start, 
	to_timestamp( st_m( st_endpoint(tr.traj) ) ) ts_end,
	ST_IsValidTrajectory(tr.traj) as valid_flag
from demo.t_taxi_trajectory tr
where tr.tid  in (28)
	order by tr.tid, tr.dt
limit 3
;
 tid |     dt     | point_count | length_kms |        ts_start        |         ts_end         | valid_flag
-----+------------+-------------+------------+------------------------+------------------------+------------
  28 | 2008-02-02 |         387 |        190 | 2008-02-02 13:33:16+08 | 2008-02-02 23:57:24+08 | t
  28 | 2008-02-03 |         885 |        331 | 2008-02-03 00:00:45+08 | 2008-02-03 23:58:18+08 | t
  28 | 2008-02-04 |         930 |        353 | 2008-02-04 00:00:09+08 | 2008-02-04 23:59:46+08 | t
(3 rows)

轨迹还原

Trajectory模型保存的数据，用ST_DumpPoints进行还原非常容易。

但是，原始GPS信息中的其他属性（速度，方向，等）无法保存和还原。

-- 还原轨迹点
with dump_pt as
(
	select (st_dumppoints(tr.traj)).geom as pt 
	from demo.t_taxi_trajectory tr
	where tr.tid=28 
		and tr.dt = '2008-02-02' 
)
select 
	st_x(pt) as lng, st_y(pt) as lat, to_timestamp(st_m(pt)) as ts
from dump_pt 
order by ts
limit 3
;
    lng    |   lat    |           ts
-----------+----------+------------------------
 116.40818 |  40.0012 | 2008-02-02 13:33:16+08
 116.42549 | 40.00747 | 2008-02-02 13:37:18+08
 116.44973 | 40.00647 | 2008-02-02 13:41:09+08
(3 rows)

PostGIS_小白笔记_构造轨迹

背景