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

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

背景技术：

球面三角形离散格网是一种经典的全球离散格网系统。它以三角形为基本格网单元，基于正八面体或正二十面体对球面空间进行同构离散化处理，形成具有层次递归性、单元形态均匀性，无缝无叠的多尺度全球格网系统(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)输入格网编码；

(2)根据编码的剖分层次l确定控制剖分层次lc；

(3)在低于控制剖分层次的阶段，采用递归逼近方法确定格网编码在控制剖分层次上所对应的控制三角形单元v0v1v2及目标编码格元在控制三角形v0v1v2内的行列号(x，y)；

(4)当转换到控制剖分层次lc时，采用直接映射法将控制三角形v0v1v2采用等分方式直接划分成格网，而后根据控制三角形内的行列号(x，y)找到目标三角形格元；

(5)计算目标三角形格元中心点坐标，将其从直角坐标转换到经纬度坐标，输出完成转换。

进一步的，所述步骤(3)中在低于控制剖分层次的阶段时，采用递归逼近法确定格网编码在控制剖分层次上所对应的控制三角形单元v0v1v2目标编码格元在控制三角形v0v1v2内的行列号(x，y)，具体过程为：

(3.1)根据格网编码获取对应初始三角形单元(例如正八面体球面三角形离散格网的初始八个球面三角形，或正二十面体球面三角形离散格网的球面三角形)，用其顶点初始化控制三角形v0v1v2；

(3.2)初始化格网剖分计数l＝0，根据各种编码与行列号的对应关系(例如二进制sqc编码，抽出其morton码的奇数位来得到的新二进制数就是列号y，morton码的偶数位组成的二进制数加上方位码就是行号x)，将格网编码转换为行列号(x，y)；

(3.3)若当前控制三角形的层次l等于控制剖分层次lc，输出控制三角形单元v0v1v2，及剩余(x，y)，递归算法结束；否则，就进入步骤(3.4)；此处剩余(x，y)是指控制三角形内的行号和列号；

(3.4)计算控制三角形单元三条弧边的中点m0、m1、m2，通过x、y判断下一层次控制三角形在当前控制三角形中的位置，若y<2^l-l-1，则为当前控制三角形的上部子三角形，进入步骤(3.5)；若x<2y-2^l-l，则为当前控制三角形的左部子三角形，令y＝y–2^l-l-1，进入步骤(3.5)；若y<2^l-l-1，则为当前控制三角形的中部子三角形，令x＝x–2×y+2^l-l-1，y＝2^l-l-1-y-1，进入步骤(3.5)；否则，为当前控制三角形的上部子三角形，令x＝x-2^l-l，y＝y-2^l-l-1，进入步骤(3.5)；

(3.5)更新控制三角形单元的顶点坐标v0、v1、v2，剖分计数l＝l+1；返回步骤(3.3)。

进一步的，所述步骤(4)中步骤(2)通过递归逼近方法获得控制剖分层次上的控制三角形单元v0v1v2后，根据剩余的剖分次数n＝l-lc,采用直接映射方法得到目标剖分层次上的三角形单元，参照平面栅格，将控制三角形v0v1v2采用等分方式直接划分成格网，而后根据控制三角形内的行列号(x，y)找到目标三角形格元，具体流程如下：

(4.1)将v0v1、v0v2等分成2ⁿ段，在v0v1上分别取第y和y+1个等分点e1、e2，在v0v2上分别取第y，y+1个等分点e3、e4；

(4.2)将e1e3等分成y段，e2e4等分成y+1段，而后在e1e3上分别取第个等分点d1、d2，在e2e4上分别取第个等分点d3、d4；

(4.3)若x为偶数，则目标层次上的三角形三个顶点为d1、d3、d4；若为奇数，则三个顶点为d4、d1、d2；

进一步的，基于步骤(4.3)所得的三个顶点计算三角形中心点坐标，并将其从直角坐标转换到经纬度坐标，完成转换。

有益效果：依据三角形格网的几何特征，通过递归逼近方法获得控制剖分层次上的控制三角形单元v0v1v2后，根据剩余的剖分次数n＝l-lc，采用直接映射方法得到目标剖分层次上的三角形单元，参照平面栅格，将控制三角形v0v1v采用等分方式直接划分成格网，而后根据控制三角形内的行列号(x，y)找到目标三角形格元。

附图说明

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

图2为本发明中基于递归逼近方法的格网编码向地理坐标的转换流程示意图；

图3为实施例中三角形格网的sqc编码与行列号对应关系示意图；

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

图5为实施例中控制三角形单元内的等分格网。

具体实施方式

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

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

(1)输入格网编码；

(2)根据编码的剖分层次l确定控制剖分层次lc；

(3)在低于控制剖分层次的阶段，采用递归逼近方法确定格网编码在控制剖分层次上所对应的控制三角形单元v0v1v2及目标编码格元在控制三角形v0v1v2内的行列号(x，y)；

(4)当转换到控制剖分层次lc时，采用直接映射法将控制三角形v0v1v2采用等分方式直接划分成格网，而后根据控制三角形内的行列号(x，y)找到目标三角形格元；

(5)计算目标三角形格元中心点坐标，将其从直角坐标转换到经纬度坐标，输出完成转换。

所述步骤(3)中在低于控制剖分层次的阶段时，采用递归逼近法确定格网编码在控制剖分层次上所对应的控制三角形单元v0v1v2目标编码格元在控制三角形v0v1v2内的行列号(x，y)，具体过程为：

(3.1)根据格网编码获取对应初始三角形单元(正八面体球面三角形离散格网的初始八个球面三角形，或正二十面体球面三角形离散格网的球面三角形)，用其顶点初始化控制三角形v0v1v2；

(3.2)初始化格网剖分计数l＝0，如图2所示，sqc编码与行列号的对应关系，将格网编码转换为行列号(x，y)；

(3.3)若当前控制三角形的层次l等于控制剖分层次lc，输出控制三角形单元v0v1v2，及剩余(x，y)，递归算法结束；否则，就进入步骤(3.4)；此处剩余(x，y)是指控制三角形内的行号和列号；

(3.4)计算控制三角形单元三条弧边的中点m0、m1、m2，通过x、y判断下一层次控制三角形在当前控制三角形中的位置，若y<2^l-l-1，则为当前控制三角形的上部子三角形，进入步骤(3.5)；若x<2y-2^l-l，则为当前控制三角形的左部子三角形，令y＝y–2^l-l-1，进入步骤(3.5)；若y<2^l-l-1，则为当前控制三角形的中部子三角形，令x＝x–2×y+2^l-l-1，y＝2^l-l-1-y-1，进入步骤(3.5)；否则，为当前控制三角形的上部子三角形，令x＝x-2^l-l，y＝y-2^l-l-1，进入步骤(3.5)；

(3.5)更新控制三角形单元的顶点坐标v0、v1、v2，剖分计数l＝l+1；返回步骤(3.3)。

所述步骤(4)中步骤(2)通过递归逼近方法获得控制剖分层次上的控制三角形单元v0v1v2后，根据剩余的剖分次数n＝l-lc,采用直接映射方法得到目标剖分层次上的三角形单元，参照平面栅格，将控制三角形v0v1v2采用等分方式直接划分成格网，如图4所示，而后根据控制三角形内的行列号(x，y)找到目标三角形格元，具体流程如下：

(4.1)将v0v1、v0v2等分成2ⁿ段，在v0v1上分别取第y和y+1个等分点e1、e2，在v0v2上分别取第y，y+1个等分点e3、e4；

(4.2)将e1e3等分成y段，e2e4等分成y+1段，而后在e1e3上分别取第个等分点d1、d2，在e2e4上分别取第个等分点d3、d4；

(4.3)若x为偶数，则目标层次上的三角形三个顶点为d1、d3、d4；若为奇数，则三个顶点为d4、d1、d2；

最后，基于步骤(4.3)所得的三个顶点计算三角形中心点坐标，并将其从直角坐标转换到经纬度坐标，完成转换。

实施例1：

本实施例中，以南京师范大学仙林校区地理科学学院所在的正二十面体三角形离散格网22层格元的sqc编码为例，使用本发明算法将其转换成地理坐标，该sqc编码为(110111010101101000010100100101011001111000000)。

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

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

步骤1：根据格网编码获取对应编号为2的初始三角形单元，用其顶点初始化控制三角形v0v1v2；

步骤2：初始化格网剖分计数l＝0，如图2所示，sqc编码与行列号的对应关系，将格网sqc编码转换为行列号(5342296，4135384)；

步骤3：若当前控制三角形的层次l等于控制剖分层次lc，输出控制三角形单元v0v1v2，及剩余行列号(x，y)(即控制三角形内的行列号)，递归算法结束；否则，就进入下一步；

步骤4：计算控制三角形单元三条弧边的中点m0、m1、m2，通过x、y判断下一层次控制三角形在当前控制三角形中的位置，若y<2^l-l-1，则为当前控制三角形的上部子三角形，进入步骤5；若x<2y-2^l-l，则为当前控制三角形的左部子三角形，令y＝y–2^l-l-1，进入步骤5；若y<2^l-l-1，则为当前控制三角形的中部子三角形，令x＝x–2×y+2^l-l-1，y＝2^l-l-1-y-1，进入步骤5；否则，为当前控制三角形的上部子三角形，令x＝x-2^l-l，y＝y-2^l-l-1，进入步骤5。

步骤5：更新控制三角形单元的顶点坐标v0、v1、v2，剖分计数l＝l+1；返回步骤3。

步骤3到步骤5是一个递归逼近过程，每次逼近行列号(x,y)的值如下：(5342296，4135384)、(1147992，2038232)、(1147992，989656)、(99416，465368)、(99416，203224)、(99416，72152)、(86183，58919)、(20647，26151)。

同时得到最终的控制三角形单元v0v1v2，其三个顶点的空间直角坐标分别为

v0(-2627.6295552955821，4766.8453899067490，3366.4147024485728)、

v1(-2599.5515195261419，4748.9227890212578，3413.2190438417356)、

v2(-2649.1708091961154，4724.8164347521324，3408.5192799219058)。

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

步骤1：将v0v1、v0v2等分成2¹⁵段，在v0v1上分别取第26151，26152个等分点e1、e2，在v0v2上分别取第26151，26152个等分点e3、e4。这几点的坐标为：

e1(-2605.2383667798254,4752.5728418204199,3403.7897166907010)、

e2(-2605.2375079911003,4752.5722913790933,3403.7911425603065)、

e3(-2644.8417093627831,4733.3408211258838,3400.0436767015644)、

e4(-2644.8423644033683,4733.3395342855165,3400.0449586184859)。

步骤2：将e1e3等分成26151段，e2e4等分成26152段，而后在e1e3上分别取第10323，10324个等分点d1、d2，在e2e4上分别取第10323，10324个等分点d3、d4。这几点的坐标为：

d1(-2620.8865472413713,4745.0080987270630,3402.3303557273480)、

d2(-2620.8880621629605,4745.0073642110137,3402.3302131326013)、

d3(-2620.8872040829283,4745.0068150879670,3402.3316399556224)、

d4(-2620.8856891612695,4745.0075496038899,3402.3317825502772)。

步骤3：因为x为奇数，则目标层次上的三角形三个顶点为d4、d1、d2；

最后，基于上述步骤所得的三个顶点计算三角形中心点坐标d(-2620.8872711624540,4745.0074260087431,3402.3307362719015)，并将其从直角坐标转换到经纬度坐标，完成转换，最终所得经纬度坐标为：p(118.91390283576838°e,32.114583311947690°n)。

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