阿里云TSDB时空数据库实战（二）-空间数据典型处理

5月5号，阿里云发布了最新的时空数据库，感兴趣的同学可以之前的文章《重磅！阿里云时空数据库正式免费公测》，
或者点击下面连接登陆阿里云官网进行免费试用：
https://www.aliyun.com/product/hitsdb_spatialpre

本文是阿里云时空数据库系列第二篇，阐述了空间数据处理中一些常见问题如何解决，为大家提供参考。第一篇是《空间数据入库与导出》

一、投影变换

提到空间数据就绕不开投影变换的问题，不少人疑惑为什么要投影变换，怎么有那么多投影，如何选择投影？
第一个问题：为什么要投影变换？首先从数据来源说起。由于我们生活的地球是一个椭球体，定位获取的数据采用地理坐标表达即经纬度，这实际是一种球体坐标，各类GPS设备之间输出这种数据；而很多几何计算需要基于平面坐标，如计算周长、面积、距离等；所以需要把球面坐标转换成平面坐标，即投影变换。
第二个问题：怎么有那么多投影？既有历史的因素，也有地球的不规则性及椭球体到平面变换的复杂因素。
第三个问题：如何选择投影？投影分为保距离、保面积、保角度三大类，还有一类是对距离、面积、角度做了均衡处理，小范围可以保证这三个衡量指标误差小；所以当选择哪种投影时需要明确保留哪个指标，如果是三个指标都要只能选择局部投影；局部投影推荐选择UTM或者高斯-克吕格。

投影带

对全球做投影带划分是为了方面数据交换、减少形变。UTM投影采用6度分带，从东经180度（或西经180度）开始，自西向东算起，因此1带的中央经线为-177（-180 -(-6)），而0度经线为30带和31带的分界，这两带的分界分别是-3和3度。中国国境所跨UTM带号为43-53 我国的疆域范围：最西端 北纬39度15分、东经73度33分 最北端 北纬53度33.5分 东经124度27分 最南点，处北纬3°51′，东经112°16′ 最东端 北纬47度27.5分 东经134度46.5分。

中国2000坐标系

中国2000坐标系又称CGCS2000，同WGS84坐标系有非常小的差异（主要是重力参数不同），如果对定位精度要求不高，如达到厘米级精度即可满足需求，可以之间拿WGS84坐标系当作CGCS2000使用。

地理坐标与平面坐标相互转换

坐标系转换对应的函数是ST_Transform
地理坐标转平面坐标示例代码如下：

SELECT ST_AsText(ST_Transform(ST_GeomFromText('POLYGON((-71.1776848522251 42.3902896512902,-71.1776843766326 42.3903829478009,
-71.1775844305465 42.3903826677917,-71.1775825927231 42.3902893647987,-71.177684
8522251 42.3902896512902))',4326),2249)) As wgs_geom;
wgs_geom
---------------------------
POLYGON((743238 2967416,743238 2967450, 743265 2967450,743265.625 2967416,743238 2967416));

平面转地理坐标系代码示例如下：

SELECT ST_AsText(ST_Transform(ST_GeomFromText('POLYGON((743238 2967416,743238 2967450,
743265 2967450,743265.625 2967416,743238 2967416))',2249),4326)) As wgs_geom;
 wgs_geom
---------------------------
POLYGON((-71.1776848522251 42.3902896512902,-71.1776843766326 42.3903829478009,-71.1775844305465 42.3903826677917,-71.1775825927231 42.3902893647987,-71.1776848522251 42.3902896512902));
(1 row)

二、线性投影

线性投影用于在一条线L上寻找到点P最近的点。整个过程包括两步：第一步点P向线L做投影，对应的函数是ST_LineLocatePoint，获得在L上的线性参考距离d；第二步利用线性参考距离d获得投影点，对应的函数是ST_LineInterpolatePoint，即为点P到L的投影点；示例代码如下：

SELECT ST_AsText(ST_LineInterpolatePoint(foo.the_line, ST_LineLocatePoint(foo.the_line, ST_GeomFromText('POINT(4 3)')))) FROM (SELECT ST_GeomFromText('LINESTRING(1 2, 4 5, 6 7)') As the_line) As foo;
st_astext
----------------
 POINT(3 4)
 

三、空间数据简化

通常使用线模型存储原始轨迹数据数据量都非常大，为了降低存储开销，可以采用损失很少的精度大幅度降低存储空间开销。利用ST_Simplify函数对线上的点做抽稀处理，抽稀的精度可以通过参数控制；示例代码如下

SELECT ST_Simplify('LINESTRING(250 250, 280 290, 300 230, 340 300, 360 260, 440 310, 470 360, 604 286)', 2);

四、几何数据合并与裁剪

通常使用线模型存储原始轨迹数据数据量都非常大，为了降低存储开销，可以采用损失很少的精度大幅度降低存储空间开销。利用ST_Simplify函数对线上的点做抽稀处理，抽稀的精度可以通过参数控制；示例代码如下

几何数据合并

当需要把两个相交或者相邻的几何体合并成一个就需要做几何体合并操作，示意图如下

这个是相交情况，经过合并处理得到

整个操作包含两步：
第一先把两个几何体使用ST_Collect函数组合在一起得到一个新几何对象
第二步再用ST_UnaryUnion函数做边界合并处理；
示例代码如下：

SELECT ST_AsText(ST_UnaryUnion(ST_Collect(st_geomfromtext('POLYGON((0 1, 0 3, 2 3, 2 1, 0 1))', 4326),
st_geomfromtext('POLYGON((1 0, 1 2, 3 2, 3 0, 1 0))', 4326))));

裁剪

用几何对象B去裁切几何对象A，保留几何对象A不与几何对象B相交的部分。示例代码如下：

SELECT ST_AsText(ST_Difference( ST_GeomFromText('LINESTRING(50 100, 50 200)'),ST_GeomFromText('LINESTRING(50 50, 50 150)')));
st_astext---------
LINESTRING(50 150,50 200)

用蓝色B对象裁剪A，裁剪后剩余的部分

五、空间数据咬合

在业务中经常遇到需要减少两个几何体缝隙的需求：比如道路边线同相邻的段之间不能有缝隙，车道面对象同相邻的车道隔离带不能有缝隙。处理方式是以其中一个几何体为主，对另外一个几何体做贴合（“咬合”），使得两个几何体之间的缝隙尽可能小。示例代码如下：

SELECT ST_AsText(ST_Snap(poly,line, ST_Distance(poly,line)*1.25)) AS polysnapped 
FROM (SELECT ST_GeomFromText('MULTIPOLYGON((( 26 125, 26 200, 126 200, 126 125, 26 125 ),( 51 150, 101 150, 76 175, 51 150 )),(( 151 100, 151 200, 176 175, 151 100 )))') 
As poly, ST_GeomFromText('LINESTRING (5 107, 54 84, 101 100)') As line
 ) As foo;

原数据线同多边形有很大的缝隙，经过上面的咬合操作得到如下结果

六、计算方位角

有些业务需要计算向量夹角，比如在交叉路口需要计算车辆右转或者左转的角度。注意下面的例子是在平面坐标下计算。
示例代码如下：

SELECT degrees(ST_Azimuth(ST_Point(25, 45), ST_Point(75, 100))) AS degA_B,degrees(ST_Azimuth(ST_Point(75, 100), ST_Point(25, 45))) AS degB_A;
  dega_b       |      degb_a
------------------+------------------
 42.2736890060937 | 222.273689006094

dega_b

dega_a