自定义经纬度索引(非RTree、Morton Code[z order curve]、Geohash的方式)

自定义经纬度索引(非RTree、Morton Code[z order curve]、Geohash的方式)

Custom Indexing for Latitude-Longitude Data

网络上有一些经纬度索引的文章,讨论的方法是用Geohash、RTree、Z-Order Curve(morton Code)等算法来将经纬度由二维变为一维的索引。

这篇文章介绍一种使用二维数据索引的方式来创建经纬度索引。前半部分讨论一种自定义的、使用二维数组的方式将经纬度变为一维数据,来加快索引。

文章最后实现了查找附近几公里的方法。常见的一些应用是根据某个经纬度点,计算出一定范围内的其他坐标点。这些坐标点是景区、商场等。中文网络里,有讲GeoHash、RTree的一些文章,但都不够全面,只是介绍了如何实现Geohash或Rtree,但没有介绍到最后一个终极目标,即给定一个经纬度 latitude\longitude和距离radius,如何实现查找一定范围内的所有的点。

本文来源于MSDN杂志的一篇文章,请参考:

https://msdn.microsoft.com/en-us/magazine/jj133823.aspx

 它有提供源代码下载,但我去下载的时候,提供的是错误的代码。

 

通常,我们在数据库里会有一张表保存如下的关系:

Id            latitude - longitude

0001      47.610 - 122.330

0002      36.175 - 115.136

0003      47.613 - 122.332

Id代表某个主键,可以是一个小吃店、商场等。Latitude – longitude 表示该ID所对应的坐标。

一些应用场景包括:

1)  要找出在一定经纬范围内的所有ID:

  SELECT Id

  FROM    IdInfoTable

  WHERE latitude >= 47.0 AND latitude < 48.0

      AND longitude >= -123.0 AND longitude < -122.0

2)还有一些应用是用户会根据自己的坐标点来查找,比如我在 47.003 – 122.311附近,我要找附近5公里范围的所有小吃店。一种最笨的方法是逐一遍历表里每一条数据,如果2个坐标点的距离小于5就把这个ID留下作为搜索结果输出。根据2个经纬度坐标点来计算2个经纬度之间距离的方法,可以参考:

http://www.cnblogs.com/softfair/p/distance_of_two_latitude_and_longitude_points.html

 

另外一种简单的想法是:

我们可以考虑把一对经纬映射到一个区域ID(sectorID)。

latitude – longitude              Sector ID

47.610 - 122.330                   49377

36.175 - 115.136                   45424

47.613 - 122.332                   49377

 

 

Id                        Sector ID

0001      49377
0002      45424
0003      49377

 

如果某些经纬度都映射到统一个SectorID,那这个ID,我都可以返回。设想有如下一个应用。我根据地图上红色坐标点,得到一个SectorID,假设为49377,我本地保存的信息表IdInfoTable里,也把经纬度信息映射为SectorID,那我做查询:

   SELECT  Id
   FROM    IdInfoTable
   WHERE   sector = 49377
就能立刻找到我需要的哪些ID是与图上红色点经纬度相同的一个区域的。

 

 

 

 

经纬度一种简单的划分,就是二维数组:

首先第一步,就是考虑如何将全球的经纬度划分开,比如我这里是每0.5度划分一小格。如下图,fraction = 0.5 表示间隔数。即每一度间隔0.5。

纬度从-90 到 90度。可以划分为:

[-90.0, -89.5)

[-89.5, -80.0)

[-80.0, -79.5)

….

[89.5, 90.0]

最后一个是闭合区间,让它包含90度。如下表格的第一列所示,这样纬度间隔总共有360度。

90-(-90)/0.5=360

经度从-180 到 180度,按照0.5间隔划分,可以表示为:

[-180.0, -179.5)

[-179.5, -170.0)

…..

[179.5, 180.0]

最后一个是闭合区间,让它包含180度。如下表格第一行所示。这样纬度间隔总共有720度。

180-(-180)/0.5=720

表格的行代表纬度,一行代表0.5度的跨度。

表格的列代表经度,一列代表0.5度的跨度。

表格总共有360 * 720 = 259,200个,序号从 0 到259,199。

fraction = 0.5

 

[-180.0, -179.5)

0

[-179.5, -170.0)

1

[-170.0, -169.5)

2

 

 

[179.5, 180.0]

719

[-90.0, -89.5) -> 0

0

1

2

...

 

719

[-89.5, -89.0) -> 1

720

721

722

...

 

1439

[-80.0, -79.5) -> 2

1440

1441

1442

...

 

 

 

...

 

 

 

 

 

 

 

 

 

 

 

 

[89.5, 90.0] -> 359

258,480

 

 

...

 

259,199

表1

 

上面使用了fraction=0.5来划分全球的经纬度。你也可以使用更小的间隔来划分全球的经纬度。这个会得到更多的sectorID。更甚者 纬度行使用0.5来划分,经度列用0.2来划分。(最好不要使用0.1这类计算机无法精确表示的浮点数。)这里演示只用0.5。

如果有一个纬度latitude的row索引 ri,有一个经度longitude的列col索引ci,则sectorID是:

sid = (ri * 360 * 1/fraction) + ci

例如:上表格里,Sector 1441的行索引 ri 是2,列索引 ci 是1。

 (2 * 360 * 1/0.5) + 1 = (2 * 720) + 1 = 1441.

360 * 1/fraction 这个因子决定了每一行有多少个列(经度360度分为多少份)。这个因子乘以ri就是这一行第一个sectorID。加上ci 就是往该行偏移ci个列,就是最终的SectorID。

 

下一个难题是如何找到经纬度的索引ri 和ci。笨的办法是遍历所有经纬度找到合适的索引,例如遍历一遍纬度,如上表格从0 到359遍历找到给定纬度的ri,再从经度的0到719遍历,找到给定经度的ci。遍历的多少取决于切分的粒度,切的越细,遍历区间越大。另外一种就是二分查找。

代码如下,纬度索引返回的是行ri值,经度索引返回的是列ci值。它们2个函数是相似的,区别只是在LonIndex函数中把常量180需要替换为360;-90 替换为-180;90替换为180。

static int LatIndex(double latitude, double fraction)

    {

      latitude = Math.Round(latitude, 8);

      int lastIndex = (int)(180 * (1.0 / fraction) - 1);

      double firstLat = -90.0;

      if (latitude == -90.0) return 0;

      if (latitude == 90.0) return lastIndex;

      int lo = 0;

      int hi = lastIndex;

      int mid = (lo + hi) / 2;

      int ct = 0; // To prevent infinite loop

      while (ct < 10000)

      {

        ++ct;

        double left = firstLat + (fraction * mid); // Left part interval

        left = Math.Round(left, 8);

        double right = left + fraction;         // Right part interval

        right = Math.Round(right, 8);

        if (latitude >= left && latitude < right)

          return mid;

        else if (latitude < left)

        {

          hi = mid - 1; mid = (lo + hi) / 2;

        }

        else

        {

          lo = mid + 1; mid = (lo + hi) / 2;

        }

      }

      throw new Exception("LatIndex no return for input = " + latitude);

    }

 

上面方法里经纬度用的double。Double、float在计算机里有精度问题,某些数据是表示不正确的,一般也不能用 = 来进行相等比较,例如 double a=-1, b=-1, 这时候计算机判断 a==b就是false。还有一些除不尽的数也是不能直接=比较的。但可以截取一个有效位数,然后进行比较。为避免equality-of-two-type-double 错误,2个double或float类型数据是不能比较大小的。所以这里取一个8位有效数字,然后进行比较。上面的函数的主循环里检查[a,b) 这个区间。所以需要显示的判断 90这个最后一个值。

有了上面的方法介绍,下面就可以把经纬度映射到sectorID了。

static int LatLonToSector(double latitude, double longitude,

  double fraction)

{

  int latIndex = LatIndex(latitude, fraction); // row

  int lonIndex = LonIndex(longitude, fraction);  // col

  return (latIndex * 360 * (int)(1.0 / fraction)) + lonIndex;

} 

 

Sector 面积的问题。

因为根据上面fraction分割经纬度的方法是将地球平面平均的分割为若干个小块。但这些小块的面积是不一样大小的。如下图:纬线之间是平行的,每一纬度之间的距离,例如A的距离、B的距离大约是111KM。经线之间是在赤道上距离远,靠近极点附近的距离近,例如C的距离是比D的距离要小的。所以一个小块的面积在不同纬度上是不一样的。

 

图1

 

 

 

2个坐标点之间的距离计算方法如下:参考:经纬度距离计算(给定2个经纬度,计算距离) http://www.cnblogs.com/softfair/p/distance_of_two_latitude_and_longitude_points.html

public static double Distance(double lat1, double lon1,

  double lat2, double lon2)

{

    double r = 6371.0; // approx. radius of earth in km

    double lat1Radians = (lat1 * Math.PI) / 180.0;

    double lon1Radians = (lon1 * Math.PI) / 180.0;

    double lat2Radians = (lat2 * Math.PI) / 180.0;

    double lon2Radians = (lon2 * Math.PI) / 180.0;

    double d = r * Math.Acos((Math.Cos(lat1Radians) *

    Math.Cos(lat2Radians) *

    Math.Cos(lon2Radians - lon1Radians)) +

      Math.Sin(lat1Radians) * Math.Sin(lat2Radians)));

    return d;

}

 

 

所以,如果知道一个Sector的起始经纬度(latitude, longitude)是能够计算这个sector的长、宽和面积的。

以下是把sectorID转化为起始纬度的函数:

static double SectorToLat(int sector,

    double fraction)

  {

    int divisor = 360 * (int)(1.0 / fraction);

    int row = sector / divisor;

    return -90.0 + (fraction * row);

  }

 

这个函数是LatLonToSector的一个逆向操作。例如给定sectorID=1441,fraction=0.5如上表格,

divisor = 360 * 1.0 / 0.5 = 720 表示有多少个列,是经度的间隔。1441 / 720 = 2 就是行row索引ri。-90.0 + 0.5 * 2 = -90.0 + 1.0 = -89.0 就是[-89.0, -79.5)这个区间的起始纬度。

获得起始经度的方法是类似的:

static double SectorToLon(int sector, double fraction)

{

  int divisor = 360 * (int)(1.0 / fraction);

  int col = sector % divisor;

  return -180.0 + (fraction * col);

}

 

sector % divisor 取模表示所在的列ci。

根据以上2个方法,可以计算面积:

static double Area(int sector, double fraction)

{

  double lat1 = SectorToLat(sector, fraction);

  double lon1 = SectorToLon(sector, fraction);

  double lat2 = lat1 + fraction;

  double lon2 = lon1 + fraction;

  double width = Distance(lat1, lon1, lat1, lon2);

  double height = Distance(lat1, lon1, lat2, lon1);

  return width * height;

}

 

 

临近区块(Adjacent Sectors)

 

有时候,找到当前Sector是不够的,还需要找到临近的Sector。例如在721块的时候,临近Sector是0, 1, 2, 720, 722, 1440, 1441 和 1442共八个块。左右Sector就是SectorId减或加1得到。上下的Sector就是SectorId减加一行所具有的列数总数即360 * (1 / fraction)。另外需要注意的几个特殊Sector,他们可能是四个角上的、或者是第一行的sector、或者是最后一行的sector,它们临近的sector是不全的。需要额外判断。以下是判断临近块的方法:

static bool IsLeftEdge(int sector, double fraction)

{

  int numColumns = (int)(1.0 / fraction) * 360;

  if (sector % numColumns == 0) return true;

  else return false;

}

static bool IsRightEdge(int sector, double fraction)

{

  if (IsLeftEdge(sector + 1, fraction) == true) return true;

  else return false;

}

static bool IsTopRow(int sector, double fraction)

{

  int numColumns = (int)(1.0 / fraction) * 360;

  if (sector >= 0 && sector <= numColumns - 1) return true;

  else return false;

}

static bool IsBottomRow(int sector, double fraction)

{

  int numColumns = (int)(1.0 / fraction) * 360;

  int numRows = (int)(1.0 / fraction) * 180;

  int firstValueInLastRow = numColumns * (numRows - 1);

  int lastValueInLastRow = numColumns * numRows - 1;

  if (sector >= firstValueInLastRow && sector <= lastValueInLastRow)

    return true;

  else

    return false;

}

 

下面的函数是返回临近点的函数。返回数组里,

result[0]是左上角的SectorId,

result[1]是正上方的SectorId,

result[2]是右上角的SectorId,

result[3]是正左方的SectorId,

result[4]是正右方的SectorId,

result[5]是左下角的SectorId,

result[6]是正下方的SectorId,

result[7]是右下角的SectorId,

 

static int[] AdjacentSectors(int sector, double fraction)

{

  int[] result = new int[8]; // Clockwise from upper-left

  int numCols = (int)(1.0 / fraction) * 360;

  int numRows = (int)(1.0 / fraction) * 180;

  int firstValueInLastRow = numCols * (numRows - 1);

  int lastValueInLastRow = numCols * numRows - 1;

  // General case

  if (IsLeftEdge(sector, fraction) == false &&

    IsRightEdge(sector, fraction) == false &&

    IsTopRow(sector, fraction) == false &&

    IsBottomRow(sector, fraction) == false)

  {

    result[0] = sector - numCols - 1;

    result[1] = sector - numCols;

    result[2] = sector - numCols + 1;

    result[3] = sector - 1;

    result[4] = sector + 1;

    result[5] = sector + numCols - 1;

    result[6] = sector + numCols;

    result[7] = sector + numCols + 1;

    return result;

  }

  // Deal with special cases here. See code download.

  throw new Exception("Unexpected logic path in AdjacentSectors");

}

 

有时候,仅仅查找临近的8个节点是不够的。例如:在高纬度地区。如下图所示,中心区域是红色的,它四周的8块Sector是绿色的。如果Fraction很小,导致红色区域的宽和高都比较小,例如有可能红色Secotor的宽高是2*10,则当需要查找附近10公里的时候,宽度需要向外延生,比如说左边,需要向外延生5个临近快。因为如果当前的点是红色Sector的起始点,即红色块左下角点时,左边就需要向外延生5*2公里。这样才能完全覆盖住半径10公里的范围。也就是下图中黑色线条标出的左边方向的5个Sector,每个Secotr的宽是2。黑色线条就一个10公里半径(以红色Sector左下角作为圆心)。

这样扩展后,会找到更多的点,包括哪些可能不是10公里范围内的点。那就需要再用距离公式对这些范围内的点再做精确的距离计算,来剔除掉一些额外点,保留下尽在范围的点。

 

 

 

 

        

       /// <summary>

        /// 根据经纬度、半径查找临近节点,会递归查找,找到该半径区域大小的所有Sector。

        /// </summary>

        /// <param name="latitude">纬度</param>

        /// <param name="longitude">经度</param>

        /// <param name="radius">半径</param>

        /// <param name="fraction"></param>

        /// <returns>SectorID</returns>

        public HashSet<int> FindClosedAdjacentSectors(double latitude, double longitude, double radius)

        {

            if (radius > 200.0)

            {

                throw new Exception("半径太大"); //超过200公里,国内查找实际意义不太大。

            }

 

            //这里使用List来保存所有的sector,再最后输出的时候再去重复。

            //因为下面的查找用的是深度优先递归查找,实际测试中发现,如果用HashSet会遗漏掉几个点。

            var allFoundSectorIdSet = new List<int>();

 

            var centralSectorId = MapLonLatCalc.LatLonToSector(latitude, longitude, Fraction);

            allFoundSectorIdSet.Add(centralSectorId);

            //计算当前sector的宽和高

            var w = 0.0; var h = 0.0;

            MapLonLatCalc.SectorArea(centralSectorId, Fraction, out w, out h);

 

            GetAdjacentSectorsRecursively(centralSectorId, allFoundSectorIdSet,radius+w,radius+h, 200);//给宽度+w,高度+h,冗余,不同经纬度上多找一些。以免遗漏

 

            return new HashSet<int>(allFoundSectorIdSet);

        }

在这里,我做递归查找时,将宽和高都额外给了一个距离,这个距离是中心块的宽和高。上面介绍的SectorToLat和SectorToLon得到的lat和lon都是Sector左上角的坐标位置(如下图蓝色点位置),一个Sector是一个矩形。如果我们输入的源点是蓝色点,那递归查找时,无需额外的增加宽和高。但如果我们输入的源点是红色的点的位置,在做半径10公里这类查找时,我们就需要额外的把红色点所在的Sector的宽和高给补足进去。代码里计算的距离的参考点都是以蓝色的点为基准的。

 

 

以上代码使用List<int> 来保存所有找到的SectorId,会存在重复的Sector,原来使用的HashSet保存,直接在查找时就去重复,但测试后发现一些问题。例如上图,中心Sector是0,首先找到周边8个Sector并保存,然后开始深度搜索递归,从Sector1首先开始,而Sector1它会优先查找到天蓝色线条所在的Sector,当找到Sector7的时候,Sector7又深度优先,会直接把Sector2找到保存下来,当Sector1周边所有都找遍完了,轮到Sector2来递归查找时,发现Sector2已经找过了,就跳过了。导致Sector8没有被包含进来。这就出现了错误。所以现在使用List保存所有Sector,最后再去重复。

 

递归查找:

 

/// <summary>

        /// 递归查找临近区域的SectorId

        /// </summary>

        /// <param name="sectorId">当前SectorID</param>

        /// <param name="allSectors">已经存在的SectorId,每找到一个就存入</param>

        /// <param name="accumulatedWidth">累计的宽度,这里其实是剩余的宽度</param>

        /// <param name="accumulatedHeight">累计的高度,这里其实剩余的宽度</param>

        /// <param name="recursiveTime">递归次数,避免死死循环</param>

        /// <param name="fraction"></param>

        private void GetAdjacentSectorsRecursively(int centorSectorId, List<int> allSectors,double remainWidth,double remainHeight, int recursiveTime)

        {

            if (recursiveTime <= 0)

            {

                return;

            }

 

            var result = MapLonLatCalc.AdjacentSectors(centorSectorId, Fraction);

 

            for (int i = 0; i < result.Length; i++)

            {

                var sId = result[i];

                if (sId > -1)

                {

                    allSectors.Add(sId);

                    //计算当前sector的宽和高

                    var w = 0.0; var h = 0.0;

                    MapLonLatCalc.SectorArea(sId, Fraction, out w, out h);

                    var leftWidth = remainWidth - w;//还剩余多少没有找到?如果>0,就需要递归查找,低纬度与高纬度值不一样。

                    var leftHeight = remainHeight - h; //看看剩余的高,高是同一经线上的2个维度距离。所以高是平均相等的。

 

                    //只要有剩余的宽和高,就递归去找。可多找,不可少找

                    // 这里有个问题,在极点附近lat接近90 的时候,Sector就是几个,无论怎么递归查找,leftWidth 和leftHeight有可能永远>0.

                    //所以在数据源上要做限制,不允许有极点附近的坐标。还好,目前国内或国外的查找点在极点附近基本没有什么。

                    if (leftWidth > 0 || leftHeight > 0)

                    {

                        GetAdjacentSectorsRecursively(sId, allSectors,leftWidth,leftHeight, recursiveTime - 1);

                    }

                }

            }

        }

在极点附近的时候,这里有个问题,在极点附近lat接近90 的时候,Sector就是几个,无论怎么递归查找,leftWidth 和leftHeight有可能永远>0.

所以在数据源上要做限制,不允许有极点附近的坐标。还好,目前国内或国外的查找点在极点附近基本没有什么。

 

测试程序:

 随机生成一批坐标点,纬度范围控制在 -85 到85之间。目前地球上85度纬度以上的地区基本无有效信息可查。微软Bing地图也在这个范围内。

static void InitLatLon()
        {
            Random r = new Random(0);

            //initialDataFile="..\\..\\UserIDLatLon.txt";
            FileStream ofs = new FileStream(initialDataFile, FileMode.Create);
            StreamWriter sw = new StreamWriter(ofs);
            for (int i = 0; i < 1000000; ++i)
            {
                double lat = (85.0 - (-85.0)) * r.NextDouble() + (-85.0);
                double lon = (180.0 - (-180.0)) * r.NextDouble() + (-180.0);
                sw.WriteLine(i.ToString().PadLeft(6, '0') + "," +
                  lat.ToString("F8") + "," + lon.ToString("F8"));
            }
            sw.Close(); ofs.Close();
        }

 

 

MyCusterInfo 是模拟现实中的一个实际业务中可能需要的一些地理信息。它除了ID、经纬度坐标外,可能还包含其他信息,这里用ExtraInfo表示。

    

public class MyCusterInfo

    {

        public int Id

        {

            get;

            set;

        }

 

 

        public double LAT

        {

            get;

            set;

        }

        public double LON

        {

            get;

            set;

        }

 

        /// <summary>

        /// 额外信息

        /// </summary>

        public object ExtraInfo

        {

            get;

            set;

        }

    }

从文件中把经纬度信息读取出来放在内存中:

                Dictionary<int, MyCusterInfo> myLatLonCusterInfo = new Dictionary<int, MyCusterInfo>(1000000);
//initialDataFile = "..\\..\\UserIDLatLon.txt";
                FileStream ifs = new FileStream(initialDataFile, FileMode.Open);
                StreamReader sr = new StreamReader(ifs);
                string line = "";
                string[] tokens = null;
                while ((line = sr.ReadLine()) != null)
                {
                    tokens = line.Split(',');

                    int id = int.Parse(tokens[0]);
                    double lat = double.Parse(tokens[1]);
                    double lon = double.Parse(tokens[2]);

                    bool vld = false;
                    double vlat = -1.0, vlon = -1.0;
                    vld = findObjects.ConvertToValidLatLon(lat, lon, out vlat, out vlon);

                    if (vld)
                    {
                        var obj=new MyCusterInfo { Id = id, LAT = vlat, LON = vlon };
                        myLatLonCusterInfo[id]=obj;
                    }
                }

 

 完整代码,请从这里获取: http://download.csdn.net/detail/soft_fair/8680673

 

经过测试程序,发现使用笨办法 和使用这个自定义经纬度索引的方法,返回的结果数是一致的。自定义索引的方法更快、更高效。目前已经应用于互联网生产环境。        

posted @ 2015-05-09 19:43 softfair 阅读(...) 评论(...) 编辑 收藏