阅读 288

PostGIS_小白笔记_构造轨迹

背景

原始的轨迹数据是一个个单独的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)
复制代码