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

一、综述

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

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

二、数据和方法

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

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

三、GeoHashUtils

    public static class Geohash{#region Direction enumpublic enum Direction{Top = 0,Right = 1,Bottom = 2,Left = 3}#endregionprivate 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;}elselon[1] = mid;}else{mid = (lat[0] + lat[1]) / 2;if (latitude > mid){ch |= Bits[bit];lat[0] = mid;}elselat[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.中关村附近道路密度