一种地理坐标向球面三角形离散格网编码的快速转换方法与流程

本发明涉及地理信息技术，具体涉及一种地理坐标向球面三角形离散格网编码的快速转换方法。

背景技术：

球面三角形离散格网是一种经典的全球离散格网系统。它以三角形为基本格网单元，基于正八面体或正二十面体对球面空间进行同构离散化处理，形成具有层次递归性、单元形态均匀性，无缝无叠的多尺度全球格网系统(tong，2013)。经过众多学者长期的研究，球面三角形离散格网不仅成为了全球多分辨率空间数据建模的有效工具，还被广泛应用于空间数据索引与可视化(goodchild，1992；feket，1990)、地图综合模型(dutton,1996)、全球导航模型(lee，2000)、全球环境监测(dutton，1999)、全球气候模拟(sadourny，1968；heikes，1995；1996；thuburn，1997)等领域。

编码作为球面三角形离散格网单元的唯一数字标识，既表达了空间位置信息，又蕴含了空间尺度信息，为处理全球多分辨率海量空间数据提供了可能，是格网系统的重要组成部分(dutton，1989，2000；fekete，1990；goodchild，1992；otoo，1993；lee，2000；bartholdi，2001)。然而，现有的大部分空间数据是以地理坐标为主，许多传统空间信息系统也是以地理坐标为基础。因此，如何准确、高效的实现球面三角形离散格网编码与地理坐标的双向转换是球面三角形离散格网研究中的核心问题之一(dutton，1999)。

由于球面是一个各向异性的二维紧致流形空间，在拓扑性质上与平面空间不同胚，属于不可展曲面(kimerling，1999)。因此，球面空间无法与平面空间一样用相同的格网连续铺砌，难以借助格网原点与目标点的距离和角度关系直接推算出目标点对应的格网编码，通常采用递归逼近的方式逐层建立地理坐标和格网编码的对应关系。在此方面，相关学者基于正八面体球面三角形格网，先后提出了zot投影(dutton，1991)、etp投影(goodchild，1992)、sqc(otoo，1993)、rca(zhao，2006)、三向互化(tong，2006)、eatp投影(sun，2007)等方法。

这些算法充分利用了基于正八面体的球面三角形格网与地理坐标具有简易转换的特性(sahr，2003)，从不同角度简化了递归逼近时的算法实现难度，提高了双向转换的效率。然而，这些算法的本质仍是基于递归逼近的方法逐层建立地理坐标与格网编码的映射关系，算法时间复杂度与格网剖分层次呈线性关系。随着格网剖分层次的提高，算法效率明显下降，不利于高分辨率空间数据的集成。同时，这些改进方法仅适用于格网结构简单但几何性质较差的正八面体球面三角形格网，难以应用在格网结构复杂但几何性质良好的正二十面体球面三角形格网上，不具有通用性。

通过分析发现，尽管现有算法充分考虑到球面空间的特殊性对球面三角形离散格网系统中格网单元编码与地理坐标双向转换的影响，但却忽略了这一特殊性质在格网剖分层次提高过程中的变化规律。球面三角形离散格网的几何性质并不是固定不变的。球面三角形离散格网在最初的剖分层次中格网单元几何性质确实差异较大，但随着剖分层次的提高，则几何性质趋于相同(white，1998)，局部区域的几何性质已近似于平面。若能根据格网几何性质的变化规律，寻找合适的控制剖分层次作为阈值，在低于该阈值的剖分层次中采用可以确保精度但效率较低的递归逼近转换方法，在高于该阈值的剖分层次采用类似平面格网的直接映射方法，则有可能构建一种兼顾精度与效率的新型混合式转换方法。

技术实现要素：

发明目的：本发明的目的在于解决现有技术中存在的不足，提供一种地理坐标向球面三角形离散格网编码的快速转换方法，其转换过程主要包括两个重要部分，控制剖分层次前的递归逼近和控制剖分层次后的直接映射。

技术方案：本发明的一种地理坐标向球面三角形离散格网编码的快速转换方法，包括以下步骤：

(1)输入经纬度坐标点p以及点p所在地图的比例尺；

(2)根据地图比例尺确定格网编码的目标剖分层次l与控制剖分层次lc；

(3)在低于控制剖分层次lc时，为保证转换精度，采用递归逼近法得到点p在控制剖分层次lc上所对应的格网行号x列号y和对应的控制三角形单元v0v1v2；

(4)当转换到控制剖分层次lc时，球面格网的曲率已趋于平面，故采用类似平面格网的直接映射方法在控制三角形v0v1v2内得到目标剖分层次上对应的最终格网行列号x和y；

(5)将最终行列号转换为格网编码，并作为结果输出。

进一步的，所述步骤(3)的具体流程如下：

(3.1)将p转换成空间直角坐标，初始化格网编码的行号x、列号y和剖分层次l：x＝0，y＝0，l＝0；

(3.2)通过点p与各初始三角形单元边界的关系确定其所在的初始三角形单元(例如：正八面体球面三角形离散格网的初始八个球面三角形，或正二十面体球面三角形离散格网的球面三角形)，并用初始三角形单元的顶点与方向初始化控制三角形顶点v0、v1、v2和方向d；

(3.3)若当前控制三角形的剖分层次l等于控制剖分层次lc，则当前行列号即为点p在控制剖分层次上对应的格网行列号，同时得到对应的控制三角形v0v1v2，递归算法结束，进入直接映射；否则，进入步骤(3.4)；

(3.4)以当前控制三角形边界v0v1的中点m0、v0v2的中点m2以及球心o构成平面om0m2，求其三维空间中的平面方程；由于平面om0m2经过球心o，则可得平面om0m2的法向量根据法方程直接给出平面om0m2的点法式方程xn×x+yn×y+zn×z＝0；将点p转换为空间直角坐标后，带入平面方程，若结果大于0，则在平面上方，下一层次的控制三角形为当前控制三角形的上部子三角形，进入步骤(3.7)；若结果小于0，则进入步骤(3.5)；

(3.5)用同样方法判断p是否在平面om1m0之上，若是，下一层次的控制三角形为当前控制三角形的左部子三角形，进入步骤(3.7)，否则进入步骤(3.6)；

(3.6)再用同样方法判断p是否在平面om2m1之上，若是，下一层次的控制三角形为当前控制三角形的右部子三角形；否则，下一层次的控制三角形为当前控制三角形的中部子三角形，进入步骤(3.7)。

(3.7)按照下表中的行列号更新规则，更新行号x和列号y，同时更新当前控制三角形块的顶点v0、v1、v2和方向d，l＝l+1，返回步骤(3.3)；：

进一步的，所述步骤(4)中，通过步骤(3)确定控制剖分层次三角形在初始三角形单元内对应的格网行列号x、y后，在余下的剖分层次n＝l-lc间采用直接映射方法获得待求点p所对应的三角形格网行列号：在控制剖分层次中确定的控制三角形单元v0v1v2内，待求点p所处的行列号tempx、tempy可由p与该控制三角形v0v1v2的相对位置求得，具体流程如下：

(4.1)在控制三角形单元v0v1v2内，过v0作v1v2的垂线交v1v2于m1，过待求点p且与v1v2平行的直线交v0m1于n1，计算列方向上的三角形格元高度进而获得控制三角形内的列号

(4.2)计算m2v2方向上的三角形格元高度从而待求点p相对于直线v0v1的偏移数结合点p所在列号tempy，点p所在的三角形格元就定位到两个格元之内，此时tempx＝2×l1(左边三角形)，或者tempx＝2×l1+1(右边三角形)；

(4.3)计算m3v1方向上的三角形格元高度从而待求点p相对于点v1的偏移数结合点p所在列号tempy，点p所在的三角形格元就定位到两个格元之内，此时tempx＝2×(l2-2ⁿ+tempy)+1(左边三角形)，或者tempx＝2×(l2-2ⁿ+tempy)+2(右边三角形)；

(4.4)步骤(4.2)和步骤(4.3)各确定两个三角形格元，它们的交集即为待求点p所在的三角形格元；先令tempx＝2×l1，若tempx≠2×(l2-2ⁿ+tempy)+1且tempx≠2×(l2-2ⁿ+tempy)+2，则令tempx＝2×l1+1；

(4.5)若控制三角形v0v1v2的方向与初始三角形方向相反，则令tempy＝2ⁿ-tempy-1，tempx＝2×tempy+tempx+1；

(4.6)finy＝y×2ⁿ+tempy，行列号finx，finy就为最终确定待求点p在初始格网中所对应的格网行列号。

有益效果：本发明依据三角形格网的几何特征，建立了在控制三角形v0v1v2内待求点p的相对位置与其对应格网编码的直接映射关系，该映射关系通过点p与控制三角形v0v1v2的相对位置关系直接定位到点p所属的目标层次三角形格元，并借助三个不同方向上的偏移量求得该格元的编码。

本发明利用球面三角形离散格网的几何性质的变化规律，在保证转换精度的前提下，能够明显提高球面三角形离散格网编码与地理坐标的转换效率，且算法时间复杂度随着剖分层次的提高，并不呈线性增长，相较于传统算法具有明显优势。

附图说明

图1为本发明的整体流程图；

图2为本发明中格网行列号的转换流程示意图；

图3为实施例中下一层次控制三角形位置判断示意图；

图4为实施例中子三角形位置分布；

图5为实施例中步骤(4.1)示意图；

图6为实施例中步骤(4.2)示意图；

图7为实施例中步骤(4.3)示意图；

图8为本发明中三角形格网的sqc编码与行列号对应关系示意图。

具体实施方式

下面对本发明技术方案进行详细说明，但是本发明的保护范围不局限于所述实施例。

如图1所示，本发明的一种地理坐标向球面三角形离散格网编码的快速转换方法，包括以下步骤：

(1)输入经纬度坐标点p以及点p所在地图的比例尺；

(2)根据地图比例尺确定格网编码的目标剖分层次l与控制剖分层次lc；

(3)在低于控制剖分层次lc时，为保证转换精度，采用递归逼近法得到点p在控制剖分层次lc上所对应的格网行号x列号y和对应的控制三角形单元v0v1v2；

(4)当转换到控制剖分层次lc时，球面格网的曲率已趋于平面，故采用类似平面格网的直接映射方法在控制三角形v0v1v2内得到目标剖分层次上对应的最终格网行列号x和y；

(5)将最终行列号转换为格网编码，并作为结果输出。

如图2所示，所述步骤(3)的具体流程如下：

(3.1)将p转换成空间直角坐标，初始化格网编码的行号x、列号y和剖分层次l：x＝0，y＝0，l＝0；

(3.2)通过点p与各初始三角形单元边界的关系确定其所在的初始三角形单元(例如：正八面体球面三角形离散格网的初始八个球面三角形，或正二十面体球面三角形离散格网的球面三角形)，并用初始三角形单元的顶点与方向初始化控制三角形顶点v0、v1、v2和方向d；

(3.3)若当前控制三角形的剖分层次l等于控制剖分层次lc，则当前行列号即为点p在控制剖分层次上对应的格网行列号，同时得到对应的控制三角形v0v1v2，递归算法结束，进入直接映射；否则，进入步骤(3.4)；

(3.4)如图3所示，以当前控制三角形边界v0v1的中点m0、v0v2的中点m2以及球心o构成平面om0m2，求其三维空间中的平面方程；由于平面om0m2经过球心o，则可得平面om0m2的法向量根据法方程直接给出平面om0m2的点法式方程xn×x+yn×y+zn×z＝0；将点p转换为空间直角坐标后，带入平面方程，若结果大于0，则在平面上方，下一层次的控制三角形为当前控制三角形的上部子三角形(如图4所示)，进入步骤(3.7)；若结果小于0，则进入步骤(3.5)；

(3.5)用同样方法判断p是否在平面om1m0之上，若是，下一层次的控制三角形为当前控制三角形的左部子三角形，进入步骤(3.7)，否则进入步骤(3.6)；

(3.6)再用同样方法判断p是否在平面om2m1之上，若是，下一层次的控制三角形为当前控制三角形的右部子三角形；否则，下一层次的控制三角形为当前控制三角形的中部子三角形，进入步骤(3.7)。

(3.7)按照下表中的行列号更新规则，更新行号x和列号y，同时更新当前控制三角形块的顶点v0、v1、v2和方向d，l＝l+1，返回步骤(3.3)；：

所述步骤(4)中，通过步骤(3)确定控制剖分层次三角形在初始三角形单元内对应的格网行列号x、y后，在余下的剖分层次n＝l-lc间采用直接映射方法获得待求点p所对应的三角形格网行列号：在控制剖分层次中确定的控制三角形单元v0v1v2内，待求点p所处的行列号tempx、tempy可由p与该控制三角形v0v1v2的相对位置求得，具体流程如下：

(4.1)如图5所示，在控制三角形单元v0v1v2内，过v0作v1v2的垂线交v1v2于m1，过待求点p且与v1v2平行的直线交v0m1于n1，计算列方向上的三角形格元高度进而获得控制三角形内的列号

(4.2)如图6所示，同理，计算m2v2方向上的三角形格元高度从而待求点p相对于直线v0v1的偏移数结合点p所在列号tempy，点p所在的三角形格元就为两个粗线格元之一，此时tempx＝2×l1(左边三角形)，或者tempx＝2×l1+1(右边三角形)；

(4.3)如图7所示，同理，计算m3v1方向上的三角形格元高度从而待求点p相对于点v1的偏移数结合点p所在列号tempy，点p所在的三角形格元就为两个粗线格元之一，此时tempx＝2×(l2-2ⁿ+tempy)+1(左边三角形)，或者tempx＝2×(l2-2ⁿ+tempy)+2(右边三角形)；

(4.4)步骤(4.2)和步骤(4.3)各确定两个三角形格元，它们的交集即为待求点p所在的三角形格元；先令tempx＝2×l1，若tempx≠2×(l2-2ⁿ+tempy)+1且tempx≠2×(l2-2ⁿ+tempy)+2，则令tempx＝2×l1+1；

(4.5)若控制三角形v0v1v2的方向与初始三角形方向相反，则令tempy＝2n-tempy-1，tempx＝2×tempy+tempx+1；

(4·6)finx＝lx/2」×2ⁿ⁺¹+tempx，finy＝y×2ⁿ+tempy，行列号finx，finy就为最终确定待求点p在初始格网中所对应的格网行列号。

最后，根据8所示的sqc编码与行列号的对应关系，将格网行列号转换为sq编码，输出结果。

实施例1：

本实施例中，以南京师范大学仙林校区地理科学学院的经纬度坐标p(118.9139003°e,32.1145815°n)转换成1:1万比例尺地图(目标剖分层次为22，控制剖分层次为7)对应的正二十面体三角形离散格网sqc编码为例。

本实施例的转换过程分为两部分，控制剖分层次前的递归逼近和控制剖分层次后的直接映射。

一、递归逼近过程步骤如下：

步骤1：将p(118.9139003°e,32.1145815°n)转换成空间直角坐标p(-2620.8871131827618,4745.0076361872807,3402.3305648447690)，初始化格网编码的行、列号和剖分层次：x＝0，y＝0，l＝0；

步骤2：通过点p与各初始三角形单元边界的关系确定其所在的初始三角形单元编号为2，并用编号为2的初始三角形单元的顶点与方向初始化控制三角形顶点v0、v1、v2和方向d；

步骤3：若当前控制三角形的剖分层次l等于控制剖分层次lc，则当前行列号即为点p在控制剖分层次上对应的格网行列号，同时得到对应的控制三角形v0v1v2，递归算法结束，进入直接映射；否则，进入下一步；

步骤4：如图2所示，以当前控制三角形边界v0v1的中点m0、v0v2的中点m2以及球心o构成平面om0m2，求其三维空间中的平面方程。由于平面om0m2经过球心o，则可得平面om0m2的法向量根据法方程可以直接给出平面om0m2的点法式方程xn×x+yn×y+zn×z＝0。将点p转换为空间直角坐标后，带入平面方程。若结果大于0，则在平面上方，下一层次的控制三角形为当前控制三角形的上部子三角形如图3所示，进入步骤7；若结果小于0，则进入下一步；

步骤5：用同样方法判断p是否在平面om1m0之上，若是，下一层次的控制三角形为当前控制三角形的左部子三角形，进入步骤7，否则进入下一步；

步骤6：再用同样方法判断p是否在平面om2m1之上，若是，下一层次的控制三角形为当前控制三角形的右部子三角形；否则，下一层次的控制三角形为当前控制三角形的中部子三角形。进入下一步。

步骤7：按照下表中的规则，更新x和y，同时更新当前控制三角形块的顶点v0、v1、v2和方向d，l＝l+1，返回步骤3。

表1行列号更新规则

步骤3到步骤7是一个递归逼近过程，每次逼近行列号(x，y)的值如下：(0，0)、(2，1)、(4，3)、(10，7)、(20,15)、(40,31)、(81,63)、(163,126)。

同时得到最终的控制三角形单元v0v1v2，其三个顶点的空间直角坐标分别为v0(-2627.6295552955821,4766.8453899067490,3366.4147024485728)、v1(-2599.5515195261419,4748.9227890212578,3413.2190438417356)、v2(-2649.1708091961154,4724.8164347521324,3408.5192799219058)，并且该三角形方向为下。

二、直接映射过程步骤如下：

步骤1：如图4所示，在控制三角形单元v0v1v2内，过v0作v1v2的垂线交v1v2于m1，过待求点p且与v1v2平行的直线交v0m1于n1。计算列方向上的三角形格元高度从而获得控制三角形内的列号

步骤2：同理，如图5所示，计算m2v2方向上的三角形格元高度从而待求点p相对于直线v0v1的偏移数

步骤3：同理，如图6所示，计算m3v1方向上的三角形格元高度从而待求点p相对于点v1的偏移数

步骤4：先令tempx＝2×l1＝10646，因为tempx≠2×(l2-2ⁿ+tempy)+1且tempx≠2×(l2-2ⁿ+tempy)+2，则令tempx＝2×l1+1＝20647。

步骤5：因为控制三角形v0v1v2的方向与初始三角形方向相反，令tempy＝2ⁿ-tempy-1＝6616，tempx＝2×tempy+tempx+1＝33880。

步骤finy＝y×2ⁿ+tempy＝4135384，行列号(5342296，4135384)就为最终确定待求点p在初始格网中所对应的格网行列号。

最后，如图8所示的sqc编码与行列号的对应关系，将格网行列号(5342296，4135384)转换为二进制的sqc格网编码(110111010101101000010100100101011001111000000)，输出结果。

本发明利用了球面三角形离散格网的几何性质的变化规律，这一规律表明球面三角形离散格网的几何性质是不断变化的，随着剖分层次的提高不断趋于近似。在低剖分层次中，由于格网单元几何性质差异较大，且格网曲率与平面差异明显，因此采用现有的层次递归逼近转换方法，以确保编码与解码的正确性。在高剖分层次中，由于格网单元几何性质近似相同，格网曲率近乎于平面，因此采用这种类似平面格网的直接映射方法，在确保转换精度的前提下，提高转换的效率。