程式隨筆
0%

快速搜尋可視範圍的座標點。 透過四元樹演算法

Posted on Edited on Disqus:

接到一個任務，要把全國的公車資料座標訊息爬回來，並且做整理歸納。研究了一下交通部的公共運輸整合資訊平台的內容，資訊相當完整，API也寫得很彈性。

唯獨公車站牌的部分，是以公車路線為出發點思考，每條路線雖然會經過同一個站牌，站牌的ID都會視為不同，不同ID的站牌，經緯度可能相同。

例如 :
富陽街口的公車站牌資料

[![](https://1.bp.blogspot.com/-rsJGRd8q0Yk/V6LuwauETRI/AAAAAAAAH2I/MXCywQPawicA3nbmx9j-D6608XWt5RxXQCLcB/s640/13765731_1287175481295263_7952967575366659415_o.jpg)](https://1.bp.blogspot.com/-rsJGRd8q0Yk/V6LuwauETRI/AAAAAAAAH2I/MXCywQPawicA3nbmx9j-D6608XWt5RxXQCLcB/s1600/13765731_1287175481295263_7952967575366659415_o.jpg)
富陽街口去程回程的公車，透過GoogleMap算大概有8支站牌，但透過API去撈會發現有17支ID皆為不同的資料。

任務目標 :
將所有站牌資料爬回來，相同的站名距離20公尺內的站牌要能自動聚合在一起，成為一支站牌。

原本思考是全部爬回來後用站名做整合，但發現一個問題是宜蘭與桃園都有相同站名的站牌”吊橋頭”，而距離相差十萬八千里，所以站名不可行。會加上20公尺的條件是因為像是台北火車站他的腹地很大，公車站牌可能很多相距很遠， 或是有些捷運站出口在兩條路口的轉角都有公車站牌，且站名相同，如果單純把這些點整合在一起放回地圖上，會發現公車站牌少了許多，失真。

歸納以上的條件，我們會有一大批公車站牌的資料(118728筆資料)，如果把每個站牌用同名稱的逐筆下去比對各自的距離，顯然掃到天荒地老也跑不完，顯然浪費效能且不切實際。

看了一篇文章後得到非常大的啟發，該方法是透過四元樹演算法將大批座標快速塞選出範圍內的點，並且做出聚合。
參考  :  How To Efficiently Display Large Amounts of Data on iOS Maps
四元樹演算法 : Wiki

節錄Wiki對四元樹的解說:「假如四元樹區塊被用來表達一組點資料(諸如一組城市的經緯度)，區塊就進行次分割直到每個葉節點包含最多一個單點。」

用圖片來表示就像這樣

[![](https://4.bp.blogspot.com/-waUqxqrmcf8/V6L7oy8YO-I/AAAAAAAAH2Y/kMcCy8J_kGw45_XPYsvcsqONAuQa767qgCLcB/s320/quadTreeInsert.gif)](https://4.bp.blogspot.com/-waUqxqrmcf8/V6L7oy8YO-I/AAAAAAAAH2Y/kMcCy8J_kGw45_XPYsvcsqONAuQa767qgCLcB/s1600/quadTreeInsert.gif)
**Source  :  https://robots.thoughtbot.com/how-to-handle-large-amounts-of-data-on-maps**
四元樹包含幾個概念 1.每個節點最多只能包含一個座標 2.當節點的座標滿了，會分裂出四個象限

根據以上概念寫出以下Class

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
/// <summary>
    /// 四元樹節點
    /// </summary>
    public class QuadTreeNode
    {
        //四個象限
        public QuadTreeNode NorthWest { get; set; }
        public QuadTreeNode NorthEast { get; set; }
        public QuadTreeNode SouthWest { get; set; }
        public QuadTreeNode SouthEast { get; set; }

        /// <summary>
        /// 邊界
        /// </summary>
        public BoundingBoxDTO BoundingBox { get; set; }

        /// <summary>
        /// 節點數量
        /// </summary>
        public int BucketCapacity { get; set; }

        /// <summary>
        /// 座標資料
        /// </summary>
        public List<Data> Points { get; set; }

        public QuadTreeNode(BoundingBoxDTO boundingBox , int capacity)
        {
            BoundingBox = boundingBox;
            BucketCapacity = capacity;
            Points = new List<BusStopAPIResult>();
        }
    }

那座標插入四元樹節點時如何運作?

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
public bool Insert(QuadTreeNode node, Data point)
        {
            // 確定這個座標符合四元樹的框框範圍
            if (!CheckBoundingBoxContainData(node.BoundingBox, point))
            {
                return false;
            }

            // 如果這個座標符合框框範圍，檢查這個節點是否已經超過座標存放上限
            // 如果沒有，那就座標放入回傳True
            if ((node.Points.Count +1) <= node.BucketCapacity)
            {
                node.Points.Add(data);
                return true;
            }

            // 會走到這一步，表示這個座標在這個四元樹的框框範圍內
            // 但是能存放的座標已經滿了，判斷這個四元樹節點是否有四個象限
            // 如果沒有，進行分裂
            if (node.NorthWest == null)
            {
                Subdivide(node);
            }

            // 將座標丟到新分裂出來的四個四元樹象限，符合座標坐落的位置，存放好回傳True
            if (Insert(node.NorthWest, data)) return true;
            if (Insert(node.NorthEast, data)) return true;
            if (Insert(node.SouthWest, data)) return true;
            if (Insert(node.SouthEast, data)) return true;

            return false;
        }

四元樹分裂的程式碼

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
/// <summary>
        /// 四元樹分裂
        /// </summary>
        /// <param name="node">四元樹節點</param>
        void Subdivide(QuadTreeNode node)
        {
            var box = node.BoundingBox;

            //找出要分裂的四元樹中心點
            double LatMid = (box.LeftTop.Latitude + box.RightBottom.Latitude) / 2.0;
            double LngMid = (box.LeftTop.Longitude + box.RightBottom.Longitude) / 2.0;

            //算出西北象限的左上右下座標點
            var NorthWest = new BoundingBoxDTO
            {
                LeftTop = new CoordinateDTO
                {
                    Latitude = box.LeftTop.Latitude,
                    Longitude = box.LeftTop.Longitude
                },
                RightBottom = new CoordinateDTO
                {
                    Latitude = LatMid,
                    Longitude = LngMid
                }
            };
            node.NorthWest = new QuadTreeNode(NorthWest, node.BucketCapacity);

            //算出東北象限的左上右下座標點
            var NorthEast = new BoundingBoxDTO
            {
                LeftTop = new CoordinateDTO
                {
                    Latitude = box.LeftTop.Latitude,
                    Longitude = LngMid
                },
                RightBottom = new CoordinateDTO
                {
                    Latitude = LatMid,
                    Longitude = box.RightBottom.Longitude
                }
            };
            node.NorthEast = new QuadTreeNode(NorthEast, node.BucketCapacity);

            //算出西南象限的左上右下座標點
            var SouthWest = new BoundingBoxDTO
            {
                LeftTop = new CoordinateDTO
                {
                    Latitude = LatMid,
                    Longitude = box.LeftTop.Longitude
                },
                RightBottom = new CoordinateDTO
                {
                    Latitude = box.RightBottom.Latitude,
                    Longitude = LngMid
                }
            };
            node.SouthWest = new QuadTreeNode(SouthWest, node.BucketCapacity);

            //算出東南象限的左上右下座標點
            var SouthEast = new BoundingBoxDTO
            {
                LeftTop = new CoordinateDTO
                {
                    Latitude = LatMid,
                    Longitude = LngMid
                },
                RightBottom = new CoordinateDTO
                {
                    Latitude = box.RightBottom.Latitude,
                    Longitude = box.RightBottom.Longitude
                }
            };
            node.SouthEast = new QuadTreeNode(SouthEast, node.BucketCapacity);

        }

看到這邊可能會有個疑問，那就是雖然知道四元樹就像是遞迴的概念，將地圖分裂成很多塊，每個大區塊可能又被分裂成若干的小區塊，一層一層的往下長

[![](https://1.bp.blogspot.com/-m_rkkVB1hxw/V6MCmq-X0gI/AAAAAAAAH2o/FHOE9VcWAvkLDp8BWJvjYnWvP9-NlbGfgCLcB/s320/500px-Point_quadtree.svg.png)](https://1.bp.blogspot.com/-m_rkkVB1hxw/V6MCmq-X0gI/AAAAAAAAH2o/FHOE9VcWAvkLDp8BWJvjYnWvP9-NlbGfgCLcB/s1600/500px-Point_quadtree.svg.png)
**Source  :  https://zh.wikipedia.org/wiki/%E5%9B%9B%E5%8F%89%E6%A0%91**
但為何會讓搜尋速度加快呢? 原因是他從上層開始看是否有重疊到要搜尋的區塊，如果有就往該節點的四象限遞迴繼續往下找重疊，如果在大區塊的節點就已經沒有包含重疊區域，那他無限往下延伸的子節點也就不會被搜尋了。引用原作者的圖片就一目了然
[![](https://3.bp.blogspot.com/-FXlbFONjFzI/V6MEINzRBHI/AAAAAAAAH20/gXhavItzQ_o3GgeyMPRgBGAxhJ9U6MuVwCLcB/s400/quadTreeQuery.gif)](https://3.bp.blogspot.com/-FXlbFONjFzI/V6MEINzRBHI/AAAAAAAAH20/gXhavItzQ_o3GgeyMPRgBGAxhJ9U6MuVwCLcB/s1600/quadTreeQuery.gif)
**Source  :  https://robots.thoughtbot.com/how-to-handle-large-amounts-of-data-on-maps**

搜尋程式碼如下
如何透過座標判斷地理位置有無交集，請參考上一篇文章: 【地理位置】透過座標框出地理位置，與判斷地理位置是否重疊

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
public void SearchDataInRange(QuadTreeNode node, BoundingBoxDTO range, ref List<Data> datas)
        {
            // 確認要搜尋的區塊是否與四元樹節點所涵蓋的區塊有重疊到，如果沒有就直接跳掉
            if (!CheckIntersectionBoundingBox(node.BoundingBox, range))
            {
                return;
            }

            //如果有，將這個區塊的節點拿出來比對是否在要搜尋的範圍中
            //因為雖然範圍有重疊到，但不代表裡面的所有節點都落在重疊的範圍
            //如果落在重疊的範圍，儲存起來，之後回傳
            foreach (var item in node.Points)
            {
                // Gather points contained in range
                if (CheckBoundingBoxContainData(range, item))
                {
                    datas.Add(item);
                }
            }

            //確認這個節點是否有分裂過
            if (node.NorthWest == null)
            {
                return;
            }

            //如果有分裂成四個象限，繼續遞迴往下找
            SearchDataInRange(node.NorthWest, range, ref datas);
            SearchDataInRange(node.NorthEast, range, ref datas);
            SearchDataInRange(node.SouthWest, range, ref datas);
            SearchDataInRange(node.SouthEast, range, ref datas);
        }

說了這麼多，我們如何將上述的四元樹理論運用在我們的目標上呢?

首先我們把一堆的公車站牌資料抓回來後，依據名稱先進行第一次分類，把每個名稱相同的依據上面的方式做成四元樹節點，但這邊有個前面提到的問題，如果是同個站名但其實是不同的站點呢?

這邊我用的方法是將同名稱的公車站牌，找出兩個離最遠的站牌算距離，如果超過1公里合理懷疑這是兩個不同的站點，只是碰巧相同名稱。台灣應該沒有差了一公里遠但是是同一個站的吧……(有的話請跟我說XD)。

將這群同站名但其實不同地方的站牌，依據一公里為單位分組，這樣就能有效地將這類同名的站牌做更細部的歸類。

計算兩點座標距離的程式

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
private const double EARTH_RADIUS = 6378.137; //地球半?

        /// <summary>
        /// 計算距離(返回公尺)
        /// </summary>
        /// <param name="lat1">點1經</param>
        /// <param name="lng1">點1緯</param>
        /// <param name="lat2">點2經</param>
        /// <param name="lng2">點2緯</param>
        /// <returns></returns>
        public static double GetDistance(double lat1, double lng1, double lat2, double lng2)
        {
            double radLat1 = rad(lat1);
            double radLat2 = rad(lat2);
            double a = radLat1 - radLat2;
            double b = rad(lng1) - rad(lng2);
            double s = 2 * Math.Asin(Math.Sqrt(Math.Pow(Math.Sin(a / 2), 2) +
             Math.Cos(radLat1) * Math.Cos(radLat2) * Math.Pow(Math.Sin(b / 2), 2)));
            s = s * EARTH_RADIUS;
            s = Math.Round(s * 10000) / 10;
            return s;
        }

        private static double rad(double d)
        {
            return d * Math.PI / 180.0;
        }

上面的問題皆排除並分類後，依據同名的站牌，算出這群站牌座標的中心點，並用離中心最遠的座標距離做為範圍，依照每小格20公尺的矩形進行四元樹搜尋，在同一格內的整合成一個站牌，記錄下來

概念圖

[![](https://2.bp.blogspot.com/-C9CqLPtlx7Q/V6MKsEK1vTI/AAAAAAAAH3E/KDH_KTeXYQQnoFWCQdLL7mvIihxNMh7KgCLcB/s640/1.png)](https://2.bp.blogspot.com/-C9CqLPtlx7Q/V6MKsEK1vTI/AAAAAAAAH3E/KDH_KTeXYQQnoFWCQdLL7mvIihxNMh7KgCLcB/s1600/1.png)
1.這些點只是亂標的，但假設真的有個同名站牌坐落在這些位置，找到中心點後，藍色線為中心點與最遠座標的直線距離，依據這個距離做為矩形二分之一寬。

2.有了距離後就能求出左上跟右下角的點，透過DbGeography這個來取得地理面積的物件

3.從左上開始，橘色的每小格為長寬各20公尺的矩形，透過四元樹下去搜尋，被同一格橘色小框框到的座標收納成同一個點，記錄下來

經過實測，如果只是雙北市的公車站點，只要4分鐘就能全部歸類掃描完畢，全國的話大概4個小時能算完，以上是這次任務的小筆記

測試環境 :
CPU  :  I7  4核 (實際在跑的時候控制大概在3核左右)
RAM : 16G