GeoHash在空间道路密度计算中的应用-以mobike骑行轨迹为例

一、综述

在空间大数据可视化展示中，往往需要计算基于道路路网的通行密度，从而得出道路不同路段的权重，然后根据不同的权重值进行道路的分色分级绘制，得出基于不同业务需求的分析展示，如我们每天日常看到的实时道路拥堵情况、道路交通流量压力等线性密度图，如下图的场景：
1.道路实时交通

2.道路通行密度（此图实际是道路的叠加，效果跟笔者所描述一致，实现方式不同）

二、数据和方法

根据摩拜的一天骑行轨迹（计算结果不受数据量的约束，理论上数据量越大越好，相应的计算量就越大），然后跟北京道路的路网中心线进行匹配计算，最终计算出哪些道路是骑行密度较大的道路，能够为政府规划、道路路网优化做决策辅助。
方法主要是根据Geohash算法做道路点和骑行轨迹点的匹配，得到道路上不同点的权重，然后根据权重进行道路点的分色渲染，最后将结果以空间化的方式展示。
摩拜骑行轨迹原始空间展示如下图所示：

具体步骤如下：
1.将路网上每条道路上的点全部提取，按照roadindex（道路索引）,pointindex（道路上点的索引）,x（当前点经度）,y（当前点维度）方式存储。
2.根据GeoHash，计算每个道路点的hash值（详细方法参考第三部分）
3.遍历摩拜骑行轨迹数据点，将每个点进行geohash转换，跟路网的每个顶点geohash匹配，如果有相同的值，则路网顶点的权重+1，最后遍历完所有的轨迹，即可计算出道路上所有的订单的权重列表
4.根据权重列表进行分级，可以按照不同的标准进行分级展示

三、GeoHashUtils

C#

    public static class Geohash
        {
         #region Direction enum

    public enum Direction
    {
        Top = 0,
        Right = 1,
        Bottom = 2,
        Left = 3
    }

    #endregion

    private const string Base32 = "0123456789bcdefghjkmnpqrstuvwxyz";
    private static readonly int[] Bits = new[] { 16, 8, 4, 2, 1 };

    private static readonly string[][] Neighbors = {
                                                       new[]
                                                           {
                                                               "p0r21436x8zb9dcf5h7kjnmqesgutwvy", // Top
                                                               "bc01fg45238967deuvhjyznpkmstqrwx", // Right
                                                               "14365h7k9dcfesgujnmqp0r2twvyx8zb", // Bottom
                                                               "238967debc01fg45kmstqrwxuvhjyznp", // Left
                                                           }, new[]
                                                                  {
                                                                      "bc01fg45238967deuvhjyznpkmstqrwx", // Top
                                                                      "p0r21436x8zb9dcf5h7kjnmqesgutwvy", // Right
                                                                      "238967debc01fg45kmstqrwxuvhjyznp", // Bottom
                                                                      "14365h7k9dcfesgujnmqp0r2twvyx8zb", // Left
                                                                  }
                                                   };

    private static readonly string[][] Borders = {
                                                     new[] {"prxz", "bcfguvyz", "028b", "0145hjnp"},
                                                     new[] {"bcfguvyz", "prxz", "0145hjnp", "028b"}
                                                 };

    public static String CalculateAdjacent(String hash, Direction direction)
    {
        hash = hash.ToLower();

        char lastChr = hash[hash.Length - 1];
        int type = hash.Length % 2;
        var dir = (int)direction;
        string nHash = hash.Substring(0, hash.Length - 1);

        if (Borders[type][dir].IndexOf(lastChr) != -1)
        {
            nHash = CalculateAdjacent(nHash, (Direction)dir);
        }
        return nHash + Base32[Neighbors[type][dir].IndexOf(lastChr)];
    }

    public static void RefineInterval(ref double[] interval, int cd, int mask)
    {
        if ((cd & mask) != 0)
        {
            interval[0] = (interval[0] + interval[1]) / 2;
        }
        else
        {
            interval[1] = (interval[0] + interval[1]) / 2;
        }
    }

    public static double[] Decode(String geohash)
    {
        bool even = true;
        double[] lat = { -90.0, 90.0 };
        double[] lon = { -180.0, 180.0 };

        foreach (char c in geohash)
        {
            int cd = Base32.IndexOf(c);
            for (int j = 0; j < 5; j++)
            {
                int mask = Bits[j];
                if (even)
                {
                    RefineInterval(ref lon, cd, mask);
                }
                else
                {
                    RefineInterval(ref lat, cd, mask);
                }
                even = !even;
            }
        }

        return new[] { (lat[0] + lat[1]) / 2, (lon[0] + lon[1]) / 2 };
    }

    public static String Encode(double latitude, double longitude, int precision = 12)
    {
        bool even = true;
        int bit = 0;
        int ch = 0;
        string geohash = "";

        double[] lat = { -90.0, 90.0 };
        double[] lon = { -180.0, 180.0 };

        if (precision < 1 || precision > 20) precision = 12;

        while (geohash.Length < precision)
        {
            double mid;

            if (even)
            {
                mid = (lon[0] + lon[1]) / 2;
                if (longitude > mid)
                {
                    ch |= Bits[bit];
                    lon[0] = mid;
                }
                else
                    lon[1] = mid;
            }
            else
            {
                mid = (lat[0] + lat[1]) / 2;
                if (latitude > mid)
                {
                    ch |= Bits[bit];
                    lat[0] = mid;
                }
                else
                    lat[1] = mid;
            }

            even = !even;
            if (bit < 4)
                bit++;
            else
            {
                geohash += Base32[ch];
                bit = 0;
                ch = 0;
            }
        }
        return geohash;
    }

java语言
参考如下网址：
https://github.com/GuoZhenghao/note/blob/master/Language/Java/Utils/GeoHash.md

四、结果

因为geohash在encode的时候需要传精度参数，所以不同精度得到的计算结果可能不一致，根据mobike一天骑行轨迹得到的骑行道路密度如下：
1.北京主城区骑行道路密度，这里划分了四个等级

2.西单附近道路密度

3.中关村附近道路密度