阅读 365

PostGIS_小白笔记_抽稀算法

背景

在处理GPS数据时,记录中往往会有很多时间戳和距离都接近的数据,一方面浪费了较多的存储空间,另一方面造成所要表达的图形不光滑或不符合标准。因此要通过某种规则,在保证矢量曲线形状不变的情况下, 最大限度地减少数据点个数,这个过程称为抽稀

常用的抽稀算法,包括但不限于(这部分内容来源于网络)

  • 步长压缩算法:沿连续曲线每隔一定的步长抽取一点,其余点全部压缩掉,然后在相邻抽取点间用直线连续或曲线拟合逼近。
  • 线段过滤算法:当某一段的长度小于某一过滤值时,就以该段的中点代替该段,如同此段的两端退化到中点一样。
  • 道格拉斯-普克(Douglas-Peuker)算法:亦称为拉默-道格拉斯-普克算法、迭代适应点算法、分裂与合并算法)是将曲线近似表示为一系列点,并减少点的数量的一种算法。

PostGIS内置了部分常用算法,我们进行一下简单的测试。

使用PostGIS函数

使用到的PostGIS函数如下

轨迹抽稀

函数ST_Simplify()基于Douglas-Peucker算法进行轨迹抽稀,具体原理这里不解释。但是有一点要注意的是,当使用WGS84时,算法的容差参数单位是度

经纬度在不同地区,每度距离差是不同的,如果假定地球是完美的球体(这样假设误差不是很大)的话,纬度为 B 的地区:

  • 纬度变化一度,球面南北方向距离变化:πR/180 ........111.7km
  • 经度变化一度,球面东西方向距离变化:πR/180cosB ....111.7cosB

比如北京 B = 40、cosB = 0.766,经度变化1度,则东西方向距离变化 85.567km。

为了简便计算,我们假设中国范围内,1度约等于100公里

Douglas-Peucker 算法

抽稀后,因为拓扑形状变化,轨迹距离会有损失。我们希望的结果是,空间比例尽量小,距离比例尽量大。

测试,当容差距离分别为1000米,100米,10米的情况下,抽稀结果的对比。

with simplified_trajectory as 
(
	-- 1 degree = 100 km
	select 
		tr.tid, tr.dt,
		st_npoints(traj) as np_orgin,
		(st_length(traj::geography)/1000) ::int as km_orgin,
		-- 0.01 degree = 1000 m
		st_npoints( ST_Simplify(traj,0.01,true) ) as np_01,
		(st_length(ST_Simplify(traj,0.01,true)::geography)/1000) ::int as km_01,
		-- 0.001 degree = 100 m
		st_npoints( ST_Simplify(traj,0.001,true) ) as np_001,
		(st_length(ST_Simplify(traj,0.001,true)::geography)/1000) ::int as km_001,
		-- 0.0001 degree = 10 m
		st_npoints( ST_Simplify(traj,0.0001,true) ) as np_0001,
		(st_length(ST_Simplify(traj,0.0001,true)::geography)/1000) ::int as km_0001,
		1 as endflag
	from demo.t_taxi_trajectory tr
)
select st.tid, st.dt, st.np_orgin, st.km_orgin,
	(st.np_01 / st.np_orgin::float)::decimal(4,2) as np_01_pct,
	(st.km_01 / st.km_orgin::float)::decimal(4,2) as km_01_pct,
	(st.np_001 / st.np_orgin::float)::decimal(4,2) as np_001_pct,
	(st.km_001 / st.km_orgin::float)::decimal(4,2) as km_001_pct,
	(st.np_0001 / st.np_orgin::float)::decimal(4,2) as np_0001_pct,
	(st.km_0001 / st.km_orgin::float)::decimal(4,2) as km_0001_pct
from simplified_trajectory st
where 1=1	
	and st.dt = '2008-02-06'
	and st.tid  in (28,10012)
	order by st.tid, st.dt
;
复制代码

Visvalingam-Whyatt 算法

ST_SimplifyVW 内置的另一个抽稀算法,没有DP那么有名。使用方法一样,这里不写语句了,有兴趣看官方文档。

测试中发现,VW算法中的容差和DP不太一样,原因未知。

步长压缩算法

对原始轨迹点使用逢十抽一的步长算法,即每10个原始点中抽一个。

-- 常规做法应该从原始GPS点数据开始,但是为了演示和测试,选择把轨迹还原成GPS然后再抽稀
with dump_pt as
(
	select tid, 
		(st_dumppoints(tr.traj)).path[1] as sn, 
		(st_dumppoints(tr.traj)).geom as pt 
	from demo.t_taxi_trajectory tr
),
-- 每10个原始点划分为一个组
rank_pt as
(
	select 
		tid, to_char(to_timestamp(st_m(pt)), 'yyyy-mm-dd') as ts ,
		st_setsrid( st_makepointm( st_x(pt), st_y(pt), st_m(pt) ), 4326) pt,
		row_number() over (partition by tid, sn/10 order by sn) as rn
	from dump_pt 
)
select rpt.tid, rpt.ts,
	st_makeline( array_agg(rpt.pt order by rpt.ts )::geometry[] ) as traj
from rank_pt rpt
where rpt.rn = 1
group by rpt.tid, rpt.ts
;
复制代码

DP算法中容差阈值的影响

  tid  |     dt     | np_orgin | km_orgin | np_01_pct | km_01_pct | np_001_pct | km_001_pct | np_0001_pct | km_0001_pct
-------+------------+----------+----------+-----------+-----------+------------+------------+-------------+-------------
    28 | 2008-02-06 |      947 |      180 |      0.01 |      0.94 |       0.03 |       0.97 |        0.15 |        0.98
 10012 | 2008-02-06 |      687 |        5 |      0.00 |      0.00 |       0.00 |       0.00 |        0.14 |        0.40
(2 rows)

复制代码

从这个样本数据中看出

  • 长原始距离(947km)
    • 容差1000米时,坐标点个数减少为1%,距离减少为94%;
    • 容差100米时,坐标点个数减少为3%,距离减少为97%;
    • 容差10米时,坐标点个数减少为15%,距离减少为98%;
  • 短原始距离(5km)
    • 容差1000米时,坐标点个数减少为0%,距离减少为0% (偏差严重);
    • 容差100米时,坐标点个数减少为0%,距离减少为0% (偏差严重);
    • 容差10米时,坐标点个数减少为15%,距离减少为40%;

所以,当原始轨迹距离比较长的情况下,用100米容差的效果很好。如果原始距离很短,需要更小的容差。

轨迹抽稀效果对比

把一条样本轨迹用不同抽稀算法处理,然后比较容量和准确度。

drop table if exists demo.t_02;
create table if not exists demo.t_02 (id text, traj geometry(linestringm, 4326));


insert into demo.t_02 
select 'src', traj from demo.t_taxi_trajectory where tid=28 and dt='2008-02-02'
union all
select 'dp', ST_Simplify(traj,0.001,true)  from demo.t_taxi_trajectory where tid=28 and dt='2008-02-02'
union all
select 'vw', ST_SimplifyVW(traj,0.00001)  from demo.t_taxi_trajectory where tid=28 and dt='2008-02-02'
;

--
with dump_pt as
(
	select tid, 
		(st_dumppoints(tr.traj)).path[1] as sn, 
		(st_dumppoints(tr.traj)).geom as pt 
	from demo.t_taxi_trajectory tr
	where tr.tid in (28) 
		and tr.dt = '2008-02-02' 
),
--
rank_pt as
(
	select 
		tid, to_char(to_timestamp(st_m(pt)), 'yyyy-mm-dd') as ts ,
		st_setsrid( st_makepointm( st_x(pt), st_y(pt), st_m(pt) ), 4326) pt,
		row_number() over (partition by tid, sn/10 order by sn) as rn
	from dump_pt 
)
--
insert into demo.t_02 
select '10-1',
	st_makeline( array_agg(rpt.pt order by rpt.ts )::geometry[] ) as traj
from rank_pt rpt
where rpt.rn = 1
group by rpt.tid, rpt.ts
;
复制代码
算法 轨迹点数量 轨迹公里数
原始轨迹 387 190
DP 0.001 130 188
VW 0.00001 107 174
步长10抽1 39 140

效果对比图如下(结论仅对本条数据有效),可见DP的准确性最好,步长10抽1最差。

拓扑抽稀

Douglas-Peucker 算法

和ST_Simplify类似,函数ST_SimplifyPreserveTopology也使用DP算法进行抽稀。但是ST_SimplifyPreserveTopology可以保证抽稀后的数据是一个合法的拓扑图形。

我不是特别明白,但是我的理解

  • 不支持M值,所以不能用来抽稀基于LinestringM的轨迹数据。
  • 是Polygon对象抽稀,最好使用这个函数。(别问为什么,因为我也不明白)
-- 容差0.001的DP拓扑抽稀
select 
	st_multi(st_simplifypreservetopology(fence_polygon, 0.001)) as fence_polygon
from demo.t_area_fence where id = 510107;
复制代码

拓扑抽稀效果对比

把一条样本围栏用不同抽稀算法处理,然后比较容量和准确度。

警告 对国家领土边界进行数据抽稀,仅仅作为验证技术可行性的研究行为,不能用于其他任何场景!!!

drop table if exists demo.t_03;
create table if not exists demo.t_03 (id text, fence_polygon geometry(multipolygon, 4326));


insert into demo.t_03 
select 'src', fence_polygon from demo.t_area_fence where id = 510107
union all
select 'dp', st_multi(ST_Simplify(fence_polygon,0.001,true)) from demo.t_area_fence where id = 510107
union all
select 'dppt', st_multi(st_simplifypreservetopology(fence_polygon, 0.001)) as fence_polygon from demo.t_area_fence where id = 510107
;
复制代码
算法 围栏点数量
原始围栏 830
DP 0.001 126
DP拓扑 0.001 127

效果对比图如下(结论仅对本条数据有效),红色方框部分有细微差异,技术研究性质的GIS业务场景下可考虑使用。

警告 对国家领土边界进行数据抽稀,仅仅作为验证技术可行性的研究行为,不能用于其他任何场景!!!