【转】在 SQL Server 2005 中使用表值函数来实现空间数据库

本文说明了如何使用 C# 和表值函数将空间搜索函数（“邻近点的点”和“多边形内的点”）添加到 Microsoft SQL Server 2005。使用此库可以在不编写任何特殊代码的情况下向应用程序中添加空间搜索。此库实现了来自约翰霍普金丝大学的公共域 C# 分层三角网格 (HTM) 算法。该 C# 库使用一组标量值函数和表值函数连接到 SQL Server 2005。这些函数起空间索引的作用。

您可以下载本文档的 Microsoft Word 版本。

本页内容

简介
表值函数：主要概念
使用表值函数添加空间索引
数据集
典型查询
结束语
附录 A：参考资料
附录 B：基本 HTM 例程

简介

空间数据搜索常用于商业和科学应用程序。结合在针对天文领域构建 SkyServer (http://skyserver.sdss.org/)方面的努力，我们开发了一种空间搜索系统。SkyServer 是一个几千 GB 大小的数据库，其中收录了大约 3 亿个天体对象。天文学家需要对该数据库进行的许多查询都会涉及到空间搜索。典型的查询包括：“邻近此点的对象有哪些”、“此区域内包含哪些对象”以及“哪些区域与此区域有重叠”？

在本文中，我们向天文学家的赤经/赤纬 (ra/dec) 天球（天空）网格添加了纬度/经度 (lat/lon) 地球网格。这两种网格大致相同，但并非精确对应，传统顺序 lat-lon 与 dec-ra 相对应。这一顺序上的颠倒迫使我们必须明确坐标系。我们将格林尼治子午线-赤道大地坐标系称为 LatLon 坐标系。该库支持三种坐标系：

• 格林尼治纬度-经度，称为 LatLon。

• 天文赤经-赤纬，称为 J2000。

• 笛卡尔 (x, y, z)，称为笛卡尔。

天文学家使用弧度分作为标准的距离度量。由于一海里为一弧度分，因此距离转换非常简单。其他许多概念也极为相似。为了说明这些问题，本文将演示如何使用此空间库针对两个 USGS 数据集构建空间索引，这两个数据集是：美国城市和美国流量计。本文使用这些索引和一些空间函数提供了若干示例，来说明如何搜索邻近某一点的城市、如何查找邻近某一城市的流量计以及如何查找某个州（多边形区域）内的流量计或城市。

我们认为此方法具有通用性。可以向几乎所有应用程序中添加空间数据核心架构和空间数据函数以便进行空间查询。这些概念也适用于其他的多维索引方案。例如，这些技术可用于搜索颜色空间或其他任何低维度度量空间。

返回页首

表值函数：主要概念

关系代数的主要概念是，每个关系运算符均使用一个或多个关系，并生成一个输出关系。SQL 是对这一概念的语法修饰，使您可以定义关系（数据定义语言）和用选择-插入-更新-删除语法来操作关系。

定义自己的标量函数使您可以对关系数据库进行扩展：您可以发送邮件，可以执行命令脚本，还可以计算非标准标量和聚合值，例如 tax() 或 median()。

但是，如果您可以创建表，则可以成为关系引擎的一部分：既是关系表的生成者，也是其使用者。这就是 OLEDB 的概念，此概念允许任何数据源生成数据流。这也是 SQL Server 2000 表值函数蕴含的概念。

在 Transact-SQL 中实现表值函数很简单：

create function t_sql_tvfPoints()
returns @points table (x float, y float)
as begin
insert @points values(1,2);
insert @points values(3,4);
return;
end

如果可以完全在 Transact-SQL 中执行函数，这样做就可以了。但是，在 SQL Server 2000 中，要在 Transact-SQL 外部实现 OLEDB 数据源或表值函数确实非常困难。

而 SQL Server 2005 集成了公共语言运行库，可以容易地创建表值函数。您可以创建列表、数组或任意 IEnumerable 对象（可以对其进行 foreach 操作的任意对象），然后将其转换为表。

[SqlFunction(TableDefinition   = "x float, y float" ,
FillRowMethodName = "FillPair")]
public static IEnumerable csPoints( )
{
int[,] points = { { 1, 2 }, { 3, 4 } };
return (IEnumerable) points;
}

在 Visual Studio 中编译这段代码，然后单击 Deploy（部署）。表值函数将被安装到数据库中。

返回页首

使用表值函数添加空间索引

对于索引，人们存在许多困惑。事实上，索引非常简单：它们只是带有一些特殊属性的表：

• SQL Server 仅有一种关联（按值）索引：B 树。B 树可以具有多字段键，但最常选择的是第一个字段。

• 从概念上说，B 树索引是由 B 树键字段、基表键字段以及您要添加到索引的任何包含字段所组成的表。

• B 树索引将根据索引键（例如邮政编码或客户 ID）来排序，以便按该键进行查找或顺序扫描会很快。

• 索引通常比基表小，只包含最重要的属性，以便与检查整个表相比，在索引中进行查找所涉及的字节数小得多。通常情况下，索引非常小，以致能完全装在主内存中，从而省去更多的磁盘存取。

• 当您执行索引查找时，可能仅仅是搜索索引（基表的垂直部分），也可能是先搜索索引，再利用基表主键将符合要求的索引逐行加入基表中（书签查找）。

中心思想是，空间索引可为您提供一小部分数据。索引会告诉您查找的位置，并通常附带一些有用的搜索信息（专家将其称为包含列或覆盖列）。索引的选择性表明了初始缩减的程度（图 1 所示的粗略子集）。找到相应的子集后，将仔细检查该子集的每个成员并舍弃假正值。图 1 中的菱形框指明了该过程。好的索引只含有少量的假正值。在整篇文章中，我们都将使用图 1 中的说法（粗略子集和仔细检查）。

图 1

B 树和表值函数可以如下组合，以使您可以构建自己的空间索引来生成粗略子集：

• 创建一个函数，以生成将相关数据聚集在一起的键。例如，如果项目 A 和 B 相关，则 A 和 B 的键在 B 树键空间中应该是邻近的。

• 创建一个表值函数，以便在给出了对所需子集的说明后，返回包含所有相关值的键范围（“覆盖”）列表。

您无法始终让每个键与其所有相关项邻近，因为键是在一个维度中排序的，而相关项是在二维或更高维数的空间中与其邻近。但是可以尽量这样做。假正值与正确答案的比率衡量了您所采取的方式的好坏。

标准方法是，找到某一条空间填充曲线，并使键空间沿该曲线穿过。例如，使用标准墨卡托地图，可以将西北中的每个人分配给西北键范围，将东南中的每个人分配给东南键范围。图 2 显示了第二种顺序空间填充曲线，它横穿所有这些象限，按顺序分配键。西北-西南象限中的每个人都具有键前缀 nwsw。如果您的区域与图 2 中所示的圆圈类似，就可以知道您的键范围

   key between 'nwsw' and 'nwse'

图 2

此搜索空间占整个表的八分之一，并且含有大约 75% 的假正值（由圆圈外但位于两个方框内的区域表示）。改进不大，但传达了一种概念。更好的索引要使用更精细的单元格分区。如果单元格足够精细，则聚合区域中的假正值就会非常少。要详细查看空间填充曲线和空间分区树，您可以参阅 Hanan Samet [Samet] 的书籍。

现在，我们要定义一种空间填充曲线：分层三角网格 (HTM)，它特别适用于球面。地球是圆的，天球也是圆的，因此，这种球面系统对于地理学家和天文学家来说非常方便。我们可以对任何度量空间做一些类似的事情。空间填充曲线提供了一些键来作为空间索引的基础。那么，如果某人具有所需的区域时，我们的表值函数将为其提供一组适当的键范围供查找（图 1 中的粗略筛选）。这些键范围将覆盖带有球面三角形的区域（称为 trixel），这与图 2 中覆盖圆圈的两个方框非常相似。搜索函数只需查看这些 trixel 的键范围内的所有对象，以确定它们是否符合要求（图 1 中的仔细检查）。

我们可以用一个具体的例子进行说明，假设有一个对象表

create table Object (objID bigint primary key,
lat   float, -- Latitude
lon   float, -- Longitude
HtmID bigint) -- The HTM key

和一个距离函数 dbo.fDistanceLatLon(lat1, lon1, lat2, lon2)，该函数可计算出两点之间的海里（弧度分）数。进一步假设以下表值函数给出了位于某个 lat-lon 点的定长半径范围内的 HtmID 点的键范围列表。

define function
fHtmCoverCircleLatLon(@lat float, @lon float, @radius float)
returns @TrixelTable table(HtmIdStart bigint, HtmIdEnd bigint)

然后，以下查询会查找旧金山 (lat,lon) = (37.8,-122.4) 周围 40 海里范围内的点。

select O.ObjID, dbo.fDistanceLatLon(O.lat,O.lon, 37.8, -122.4)
from fHtmCoverCircleLatLon(37.8, -122.4, 40) as TrixelTable
join Object O
on O.HtmID between TrixelTable.HtmIdStart        -- coarse test
and TrixelTable.HtmIdEnd
where dbo.fDistanceLatLon(lat,lon,37.8, -122.4) < 40   -- careful test

现在，我们必须定义 HTM 键生成函数、距离函数和 HTM 覆盖函数。下一步我们将以两组美国地质空间数据集为例执行这些操作。如果您不相信其中包含的对象达数十亿，请访问 http://skyserver.sdss.org/ 并浏览一下该网站。该网站使用相同的代码来对几千 GB 的天文数据库进行空间查找。

本文主要讲述如何使用 SQL 表值函数和像 HTM 这样的空间填充曲线来构建空间索引。同样地，我们将 HTM 代码本身当作一种在别处 [Szalay] 存档的黑色方框，我们只关注如何使其在 SQL 应用程序内适合我们的需要。

返回页首

数据集

美国地质勘探局 (USGS) 收集并发布了美国的相关数据。图 3 显示了由 USGS 维护的 18000 台用于测量河流水流量和级别的流量计。USGS 还发布了 23000 个地点的地名及其人口的列表。

图 3

USGS 居民点（23000 个城市）

USGS 在 1993 年发布了地点名称及其某些属性的列表。USGS 网站上有一些更新的列表，但由于它们是按州来分段的，因此很难获得一个全国范围的列表。旧的列表足以能够演示空间索引。数据格式如下：

create table Place(
PlaceName   varchar(100) not null, -- City name
State       char(2)      not null, -- 2 char state code
Population  int          not null, -- Number of residents (1990)
Households  int          not null, -- Number of homes (1990)
LandArea    int          not null, -- Area in sqare KM
WaterArea   int          not null, -- Water area within land area
Lat         float        not null, -- Latitude in decimal degrees
Lon         float        not null, -- Longitude decimal degrees
HtmID       bigint       not null primary key --spatial index key
)

为了加快名称查找，我们添加了一个名称索引，但数据是按空间键聚集在一起的。邻近的对象共同位于聚集 B 树中，并因此位于相同或相邻的磁盘页上。

create index Place_Name on Place(PlaceName)

可以从 USGS 网站下载除 HtmID 数据以外的其他所有数据。可以用 SQL Server 2005 数据导入向导来导入数据（我们已在示例数据库中进行了此操作）。HtmID 字段是根据 LatLon 按以下方法计算出的：

update Place set HtmID = dbo.fHtmLatLon(lat, lon)

USGS 流量计（17000 台仪器）

USGS 自 1854 起一直在维护河流流量记录。到 2000 年 1 月 1 日，他们已累积了超过 43 万年的测量数据。大约有六千个活动的测量站处于活动状态，而有大约四千个处于联机状态。http://waterdata.usgs.gov/nwis/rt 中对这些流量计进行了详细介绍。有个 NOAA 站点以非常便利的方式显示了来自几百个最普遍的测量站的数据：http://weather.gov/rivers_tab.php。

我们的数据库只包含美国大陆的测量站（见图 3）。关岛、阿拉斯加、夏威夷、波多黎各和维京群岛也有测量站，但不包含在此数据库内。流量计测量站表为：

create table Station (
StationName   varchar(100) not null,     -- USGS Station Name
State         char(2)      not null,     -- State location
Lat           float        not null,     -- Latitude in Decimal
Lon           float        not null,     -- Longitude in Decimal
DrainageArea  float        not null,     -- Drainage Area (km2)
FirstYear     int          not null,     -- First Year operation
YearsRecorded int          not null,     -- Record years (at Y2k)
IsActive      bit          not null,     -- Was it active at Y2k?
IsRealTime    bit          not null,     -- On Internet at Y2K?
StationNumber int          not null,     -- USGS Station Number
HtmID         bigint       not null,     -- HTM spatial key
-- (based on lat/lon)
primary key(htmID, StationNumber) )

如上所述，HtmID 字段是根据 LatLon 按以下方法计算出的：

update Station set HtmID = dbo.fHtmLatLon(lat, lon)

由于一个位置有多达 18 个测量站，因此主键必须包括测量站编号以便区分它们。但是，在 B 树中，HTM 键将所有邻近的测量站聚集在一起。为了加快查找，我们添加了测量站编号和名称索引：

create index Station_Name   on Station(StationName)
create index Station_Number on Station(StationNumber)

空间索引表

现在，我们就可以创建自己的空间索引了。我们本可以将字段添加到基表，但要使存储过程对多个不同的表均有效，我们发现，只需将所有对象加入到一个空间索引中即可。您可以选择 (type,HtmID) 作为键来隔离不同类型的对象；我们选择了 (HtmID, key) 作为键，以便将所有类型（城市和流量计）的邻近对象聚集在一起。该空间索引为：

create table SpatialIndex (
HtmID   bigint   not null , -- HTM spatial key (based on lat/lon)
Lat     float    not null , -- Latitude in Decimal
Lon     float    not null , -- Longitude in Decimal
x       float    not null , -- Cartesian coordinates,
y       float    not null , -- derived from lat-lon
z       float    not null , --,
Type    char(1)  not null , -- Place (P) or gauge (G)
ObjID   bigint   not null , -- Object ID in table
primary key (HtmID, ObjID) )

此主题后面将对笛卡尔坐标进行说明。至于现在，我们只需要知道，fHtmCenterPoint(HtmID) 函数将返回该 HTM 三角形中心点的笛卡尔 (x,y,z) 单位向量。这就是该 HTM 的极限点，因为此中心被细分为无穷个小的 trixel。

SpatialIndex 表将根据 Place 和 Station 表的数据来填充，如下所示：

insert SpatialIndex
select   P.HtmID, Lat, Lon, XYZ.x, XYZ.y, XYZ.z,
'P' as type, P. HtmID as ObjID
From   Place P cross apply fHtmLatLonToXyz(P.lat, P.lon)XYZ
insert SpatialIndex
select S.HtmID, Lat, Lon, XYZ.x, XYZ.y, XYZ.z,
'S' as type, S.StationNumber as ObjID
from   Station S cross apply fHtmLatLonToXyz(S.lat, S.lon) XYZ

为了清理数据库，应执行以下命令：

DBCC INDEXDEFRAG    (spatial  , Station, 1)
DBCC INDEXDEFRAG    (spatial  , Station, Station_Name)
DBCC INDEXDEFRAG    (spatial  , Station, Station_Number)
DBCC INDEXDEFRAG    (spatial  , Place,   1)
DBCC INDEXDEFRAG    (spatial  , Place,   Place_Name)
DBCC INDEXDEFRAG    (spatial  , SpatialIndex, 1)
DBCC SHRINKDATABASE (spatial  , 1  ) - 1 percent spare space

题外话：笛卡尔坐标

您可以选择跳过此部分。此部分对于使用该库并不是必需的。HTM 代码必须依靠一种技巧来避开球面几何问题：它从球面的二维表面移到了三维空间。这就使得“在多边形内”和“在点附近”查询的检查进行得非常快。

球面上的每个 lat/lon 点都可以表示为三维空间中的一个单位向量 v = (x,y,z)。北极和南极（90° 和 -90°）的单位向量分别为 v = (0,0,1) 和 v = (0,0,-1)。Z 代表旋转轴，XZ 平面代表本初（格林尼治）子午线，其经度为 0° 或 180°。正式的定义为：

x = cos(lat)cos(lon)
y =cos(lat)sin(lon)
z = sin(lat)

这些笛卡尔坐标的使用方法如下：给定单位球面上的两点 p1=(x1,y1,z1) 和 p2 = (x2,y2,z2)，则它们的点乘积 p1&p2 = x1*x2+y1*y2+z1*z2 就是这两点所代表的单位向量之间的角度的余弦值。它是一个距离度量。图 4 显示了笛卡尔坐标如何使得对“多边形内的点”和“邻近点的点”的检查快速进行。每个 lat/lon 点均具有一个对应的 (x,y,z) 单位向量。

图 4

如果我们要查找 p1 点周围 45 海里（弧度分）范围内的点，则它与 p1 点最多成 45/60 度。这些点与 p1 的点乘积将小于 d=acos(radians(45/60)。该“邻近”检查即变为 { p2 | p2&p1 < d}，它将很快进行。在图 5 中，每个大圆圈或小圆圈都是一个平面与该圆圈的交集，如果某一点与平面法向量的点乘积小于 cos(?¨)（其中 2?¨ 是该圆圈的弧度角直径），则该点就在圆圈内部。

图 5

笛卡尔坐标还可以使得对“多边形内的点”的检查快速进行。所有的多边形都具有大圆边或小圆边。这些边沿着与球面交叉的某个平面分布。因此，可以通过与该平面垂直的单位向量 v 及其位移来定义这些边。例如，赤道就是向量 v = (0,0,1)，且位移为零。纬度 60° 由向量 v = (0,0,1) 及位移 0.5 来定义，而围绕巴尔的摩市的 60° 纬度圈由向量 v = (0.179195, -0.752798, 0.633392) 及 0.5 个位移来定义。对于地点 p2，如果 p2&v < 0.5，则该地点就位于巴尔的摩市 60° 纬度圈内。同样，可以通过计算出三个这种点乘积来确定某个点位于 HTM 三角形内部还是外部。这是 HTM 代码如此有效和快捷的主要原因之一。

我们实现了若干辅助过程来进行从 LatLon 到笛卡尔坐标的转换。

本文稍后将用到这些函数，它们存档在 API spec and Intellisense [Fekete] 中。

在此，该库的默认设置为 21 级深度的 HTM 键（第一级将球面分成八个表面区域，随后每一级将球面三角形分成四个子三角形）。从表 1 中可看出，21 级深度的 trixel 相当小。最多可以将代码修改为 31 级深度，这是因为不能超出 64 位表示法所占用的位数。

在表 1 中，每个 HTM 级别都会细分球面。对于每个级别，此表均以度数、弧度分、弧度秒和米四种单位的平方形式表示了面积。trixel 列显示了一些特征大小：默认的 21 级深度的 trixel 大约为 0.3 平方弧度秒。USGS 数据在每 12 级深度的 trixel 具有大约半个对象。

表 1

HTM 深度			面积				对象/trixel
	deg 2	arc min 2	arc sec 2	earth m 2	trixel	SDSS	USGS
球面	41253	148,510,800	534,638,880,000	5.E+14
0	5157	18,563,850	66,829,860,000	6E+13		3E+8
1	1289	4,640,963	16,707,465,000	2E+13		8E+7
2	322	1,160,241	4,176,866,250	4E+12		2E+7
3	81	290,060	1,044,216,563	1E+12		5E+6
4	20	72,515	261,054,141	2E+11		1E+6	30,000
5	5	18,129	65,263,535	6E+10		3E+5	7,500
6	1	4,532	16,315,884	2E+10	1 deg2	73242	1,875
7	3E-1	1,133	4,078,971	4E+9		18311	468
8	8E-2	283	1,019,743	1E+9		4578	117
9	2E-2	71	254,936	2E+8		1144	29
10	5E-3	18	63,734	6E+7		286	7
11	1E-3	4	15,933	2E+7		72	2
12	3E-4	1	3,983	4E+6	1 amin2	18	0.5
13	8E-5	3E-1	996	943816		4	0.1
14	2E-5	7E-2	249	235954		1
15	5E-6	2E-2	62	58989		0.3
16	1E-6	4E-3	16	14747		.
17	3E-7	1E-3	4	3687
18	8E-8	3E-4	1	922
19	2E-8	7E-5	2E-1	230	1 asec2
20	5E-9	2E-5	6E-2	58	1 km2
21	1E-9	4E-6	2E-2	14
22	3E-10	1E-6	4E-3	4
23	7E-11	3E-7	9E-4	1	1 m2
24	2E-11	7E-8	2E-4	2E-1
25	5E-12	2E-8	6E-5	6E-2
26	1E-12	4E-9	1E-5	1E-2

返回页首

典型查询

现在，您应当可以执行一些查询了。

查询 1：查找邻近点的点：查找邻近某一地点的城镇

最常见的查询是查找邻近某个特定地点或点的所有地点。考虑以下查询，“查找马里兰州巴尔的摩市周围 100 海里内的所有城镇”。通过下面的方法来得到覆盖巴尔的摩市周围 100 海里（100 弧度分）范围的 HTM 三角形

select * -- find a HTM cover 100 NM around Baltimore
from fHtmCoverCircleLatLon(39.3, -76.6, 100)

它将返回表 2 中所示的 Trixel 表。即，fHtmCoverCircleLatLon() 函数将返回“覆盖”该圆圈（在本例中，是单个 trixel）的一组 HTM 三角形。该圆圈内所有对象的 HTM 键也位于这些三角形中之一内。现在，我们需要检查所有这些三角形并舍弃假正值（图 1 中的仔细检查）。我们将按照与巴尔的摩市的距离对答案集进行排序，因此，如果我们需要找出最近的地点，只需选择 TOP 1 WHERE distance > 0 即可（我们要从中排除巴尔的摩市本身）。

declare @lat float, @lon float
select @lat = lat, @lon = lon
from Place
where Place.PlaceName = 'Baltimore'
and State = 'MD'
select ObjID, dbo.fDistanceLatLon(@lat,@lon, lat, lon) as distance
from SpatialIndex join fHtmCoverCircleLatLon(@lat, @lon, 100)
On HtmID between HtmIdStart and HtmIdEnd           -- coarse test
and type = 'P'
and dbo.fDistanceLatLon(@lat,@lon, lat, lon) < 100 -- careful test
order by distance asc

表 2. 巴尔的摩市 HTM 覆盖纬度圈

HtmIdStart	HtmIdEnd
14023068221440	14027363188735

此覆盖联接将返回 2928 行（粗略检查）；其中 1122 行在 100 航空英里以内（仔细检查）。它给出了 61% 的假正值：全部操作在 9 毫秒内完成。

由于这些是常见的任务，因此具有针对它们的标准函数：

fHtmNearbyLatLon(type, lat, lon, radius)
fHtmNearestLatLon(type, lat, lon)

这样，上述查询就变为

select ObjID, distance
from fHtmNearestLatLon('P', 39.3, -76.61)

查询 2：查找某个方框内的地点

在显示正方形的地图或窗口时，应用程序通常需要查找某个正方形视区内的所有对象。科罗拉多州几乎完全是正方形的，它的西北角点为 (41°N-109°3'W)，西南角点为 (37°N-102° 3'E)。该州的中心点为 (39°N, -105°33'E)，因此可以用以该点为中心的圆圈覆盖该正方形。

declare @radius float
set @radius = dbo.fDistanceLatLon(41,-109.55,37,-102.05)/2
select *
from Station
where StationNumber in (
select ObjID
from fHtmCoverCircleLatLon(39, -105.55, @radius) join SpatialIndex
on HtmID between HtmIdStart and HtmIdEnd
and lat between 37 and 41
and lon between -109.05 and -102.048
and type = 'S')
OPTION (FORCE ORDER)

本例在大约 46 毫秒内返回了 1030 个流量计。科罗拉多州的其他五个流量计正好在其边界上，散布于距离南纬 37°纬度圈不超过一海里的范围内。如果将南纬从 37° 调整到 36.98°，则其他这五个测量站就会出现。（GIS 系统和天文应用程序通常需要某一区域周围具有缓冲区。此 HTM 代码包含对缓冲区的支持，在实际的应用程序经常会用到缓冲区。请查看参考资料 [Szalay] 以了解具体操作方式。）此覆盖圆圈返回了 36 个三角形。与 SpatialIndex 表的联接返回了 1975 个流量计。其中包含 47% 的假正值。下一节将演示如何通过使用 HTM 区域指定多边形覆盖而非圆圈覆盖对此进行改进。

FORCE ORDER 子句比较麻烦：如果缺少该子句，查询的运行时间会长十倍，因为优化器会将空间索引作为外部表进行嵌套循环联接。如果这些表更大（包含数百万行），优化程序有可能会采取其他计划，但我们不能指望它。优化程序不可能无需提示就能针对上一部分中的所有查询选择正确的计划。

查询 3：查找某个多边形内的地点

HTM 代码允许我们将此区域指定为圆圈、矩形、凸球面或这些区域的组合。特别地，HTM 库允许使用下面的线性语法来指定区域：

circleSpec  := 'CIRCLE LATLON '      lat lon radius
|  'CIRCLE J2000 '       ra  dec radius
|  'CIRCLE [CARTESIAN ]' x y z   radius
rectSpec    := 'RECT LATLON '        { lat lon }2
|  'RECT J2000 '         { ra  dec }2
|  'RECT [CARTESIAN ]'   { x y z   }2
hullSpec    := 'CHULL LATLON '       { lon lat }3+
|  'CHULL J2000 '        { ra dec  }3+
|  'CHULL [CARTESIAN ]'  { x y z   }3+
convexSpec  := 'CONVEX ' [ 'CARTESIAN '] { x y z D }*
areaSpec    := rectSpec | circleSpec | hullSpec | convexSpec
regionSpec  := 'REGION ' {areaSpec}* | areaSpec

下面给出了区域指定示例：

• 圆圈。指定一个点和大小为 1.75 海里（弧度分）的半径。

  'CIRCLE LATLON 39.3 -76.61 100'
      'CIRCLE CARTESIAN 0.1792 -0.7528 0.6334 100'
      

• 矩形。指定用来定义最小和最大 lat,lon 的两个角点。经度坐标以折回方式确定，即 lonmin=358.0 和 lonmax=2.0，这是一个四度宽的范围。纬度必须介于北极和南极之间。矩形边是纬度或经度保持不变的直线，而非凹形和凸形的大圆边。

  'RECT LATLON 37 -109.55  41 -102.05'
      

• 凹形。指定用来定义凸球面三个或更多个点，且该凸球面的边用大圆圈将相邻的点连接起来。这些点必须位于单个半球内，否则会发生错误。这些点的顺序无关紧要。

  'CHULL LATLON 37 -109.55 41 -109.55 41 -102.051 37 -102.05'
      

• 凸形。以笛卡尔向量 (x,y,z) 及该向量单位长度的分量的形式指定任意多个（包括零个）约束。

  'CONVEX   -0.17886 -0.63204 -0.75401 0.00000
      -0.97797  0.20865 -0.00015 0.00000
      0.16409  0.57987  0.79801 0.00000
      0.94235 -0.33463  0.00000 0.00000'
      

• 区域。区域是零个或更多个圆圈、矩形、凹形和凸形区域的组合。

  'REGION CONVEX 0.7 0.7 0.0 -0.5 CIRCLE LATLON 18.2 -22.4 1.75'
      

可以这些区域描述中的任意一个应用于 fHtmCoverRegion() 例程，以返回一个 trixel 表，用来描述覆盖该区域的一组 trixel（三角形区域）。用于此科罗拉多州查询的较简单的代码为：

select S.*
from (select ObjID
from fHtmCoverRegion('RECT LATLON 37 -109.55  41 -102.05')
loop join SpatialIndex
on HtmID between HtmIdStart and HtmIdEnd
and lat between 37 and 41
and lon between -109.05 and -102.048
and type = 'S') as G
join Station S on G.objID = S.StationNumber
OPTION (FORCE ORDER)

有必要使用这种不寻常的查询格式来准确告诉优化器联接的执行顺序（以使“强制顺序”选项正确运行）。很难以这种方式修改优化器，直到表值函数具有统计信息为止，估计它们会非常耗费资源。您必须强制使它们进入内部循环联接。

此查询将返回 1030 个流量计，而有 1365 个来自覆盖范围的候选项，因此包含 25% 的假正值。请注意，矩形覆盖优于圆形覆盖，因为后者包含 61% 的假正值。对于非矩形形状的州，可以使用多边形语法，但本文仅讲述表值函数，而不讲述 HTM 算法。您可以在项目中以及项目的相关文档中查看 HTM 代码。

可以转换为凸球面进行类似的查询。

select S.*
from (select ObjID
from fHtmCoverRegion(
'CHULL LATLON 37 -109.55 41 -109.55 41 -102.05 37 -102.05')
loop join SpatialIndex
on HtmID between HtmIdStart and HtmIdEnd
and lat between 37 and 41
and lon between -109.05 and -102.048
and type = 'S') as G
join Station S on G.objID = S.StationNumber
OPTION (FORCE ORDER)

此查询将返回 1030 个流量计，而有 1193 个来自覆盖范围的候选项，因此包含 14% 的假正值。在本例中，凸球面覆盖比同等的矩形覆盖更好。

查询 4：高级主题：复杂区域

前面的示例给出了用于区域的语法，并对“邻近点的点”和“矩形内的点”搜索进行了论述。区域可能会十分复杂。它们是多个凸形区域的布尔组合。我们无法在此详细解释区域的概念，但伴随项目中的 HTM 库包含对区域进行布尔组合、简化区域、计算区域顶点和计算区域面积的逻辑，还包含许多其他特性。[Fekete]、[Gray] 和 [Szalay] 中介绍了这些概念。

为了初步理解这些概念，我们以犹他州为例。它的边界可用两个矩形来近似地定义：

declare @utahRegion varchar(max)
set @utahRegion = 'region '
+ 'rect latlon 37 -114.0475  41 -109.0475 ' -- Main part
+ 'rect latlon 41 -114.0475  42 -111.01  '  -- Ogden and Salt Lake.

现在，我们可以用以下查询来查找犹他州中的所有流量计：

select S.*
from (
select ObjID
from fHtmCoverRegion(@utahRegion)
loop join SpatialIndex
on HtmID between HtmIdStart and HtmIdEnd
and (((     lat between 37        and      41)  -- Careful test
and (lon between -114.0475 and -109.04)) -- Are we inside
or ((   lat between 41        and      42)  -- one of the two
and (lon between -114.0475 and -111.01)) -- boxes?
)
and type = 'S' ) as G
join Station S on G.objID = S.StationNumber
OPTION (FORCE ORDER)

覆盖返回了 38 个 trixel。联接返回了 775 个测量站。仔细检查找到了犹他州中的 670 个测量站，另外有怀俄名州的两个测量站正好位于接壤处（14% 的假正值）。

大多数州需要使用复杂得多的区域。例如，近似地将加利福尼亚州的边界连起来的区域为：

declare @californiaRegion varchar(max)
set @californiaRegion = 'region '
+ 'rect latlon 39    -125 '    -- Nortwest corner
+   '42    -120 '    -- Center of Lake Tahoe
+ 'chull latlon 39    -124 '    -- Pt. Arena
+ '39    -120 '    -- Lake Tahoe.
+ '35    -114.6 '  -- Start Colorado River
+ '34.3  -114.1 '  -- Lake Havasu
+ '32.74 -114.5 '  -- Yuma
+ '32.53 -117.1 '  -- San Diego
+ '33.2  -119.5 '  -- San Nicholas Is
+ '34    -120.5 '  -- San Miguel Is.
+ '34.57 -120.65 ' -- Pt. Arguelo
+ '36.3  -121.9 '  -- Pt. Sur
+ '36.6  -122.0 '  -- Monterey
+ '38    -123.03 ' -- Pt. Rayes
select stationNumber
from fHtmCoverRegion(@californiaRegion)
loop join SpatialIndex
on HtmID between HtmIdStart and HtmIdEnd
/* and <careful test> */
and type = 'S'
join Station S on objID = S.StationNumber
OPTION (FORCE ORDER)

覆盖返回了 108 个 trixel，一共包含 2132 个测量站。其中，1928 个位于加利福尼亚州内，因此假正值比率大约为 5%。但是仔细检查并非小事。

针对地点而非测量站进行的相同查询（包括仔细检查）类似于以下代码：

select *
from Place
where HtmID in
(select distinct SI.objID
from fHtmCoverRegion(@californiaRegion)
loop join SpatialIndex SI
on SI.HtmID between HtmIdStart and HtmIdEnd
and SI.type = 'P'
join place P on SI.objID = P.HtmID
cross join fHtmRegionToTable(@californiaRegion) Poly
group by SI.objID, Poly.convexID
having min(SI.x*Poly.x + SI.y*Poly.y + SI.z*Poly.z - Poly.d) >= 0
)
OPTION (FORCE ORDER)

此查询使用加利福尼亚州的凸形半空间表示法和 [Gray] 中介绍的技术来快速检查某个点是否位于加利福尼亚州凸球面内。它返回了 885 个地点，其中 7 个位于与加利福尼亚州毗邻的亚利桑那州（多边形近似于加利福尼亚州的边界）。它在 1GHz 的处理器上运行了 0.249 秒。如果不用 OPTION(FORCE ORDER) 子句，其运行速度将变慢，需要花费 247 秒。

由于此要求十分常见，而且代码极具技巧性，因此我们添加了 fHtmRegionObjects(Region,Type) 过程，用来从空间索引中返回对象 ID。由于此过程封装了前面所示的两个技巧性代码，因此这两个针对加利福尼亚州的查询变为了：

select *                 -- Get all the California River Stations
from Station
where stationNumber in   -- that are inside the region
(select ObjID
from fHtmRegionObjects(@californiaRegion,'S'))
select *                 -- Get all the California cities
from Place
where HtmID in           -- that are inside the region
(select ObjID
from fHtmRegionObjects(@californiaRegion,'P'))

针对科罗拉多州和犹他州的查询也可以使用此例程来简化。

返回页首

结束语

在此所述的 HTM 空间索引库本身是有趣而又有用的。基于球面为“多边形内的点”查询索引数据是一种便利的方法。但是，该库还作为一个示例很好地说明了如何通过添加以诸如 C#、C++、Visual Basic 或 Java 之类的语言进行实际计算的类库来扩充 SQL Server 和其他数据库系统。实现功能强大的表值函数和标量函数并将这些查询及其永久数据集成到数据库中的能力是一种非常强大的扩充机制，它将在保证对象关系数据库的基础上传递。这仅仅是个开头。在接下来的十年中，编程语言和数据库查询语言有可能获得更好的数据集成方式。这对于应用程序开发人员来说将是一件好事。

详细信息请参见：

http://msdn.microsoft.com/sql/

项目编辑：Susanne Bonney

返回页首

附录 A：参考资料

• [Gray]“There Goes the Neighborhood:Relational Algebra for Spatial Data Search”。Jim Gray、Alexander S. Szalay、Gyorgy Fekete、Wil O'Mullane、Maria A. Nieto-Santisteban、Aniruddha R. Thakar、Gerd Heber、Arnold H. Rots，MSR-TR-2004-32，2004 年 4 月

• [Szalay]“Indexing the Sphere with the Hierarchical Triangular Mesh”。Alexander S. Szalay、Jim Gray、George Fekete、Peter Z. Kunszt、Peter Kukol、Aniruddha R. Thakar，Microsoft SQL Server 2005 Samples。

• [Fekete]“SQL SERVER 2005 HTM Interface Release 4”。George Fekete、Jim Gray、Alexander S. Szalay，2005 年 5 月 15 日，Microsoft SQL Server 2005 Samples。

• [Samet1]“Applications of Spatial Data Structures:Computer Graphics, Image Processing, and GIS”。Hanan Samet，Addison-Wesley, Reading, MA, 1990。ISBN0-201-50300-0。

• [Samet2]“The Design and Analysis of Spatial Data Structures”。Hanan Samet，Addison-Wesley, Reading, MA, 1990。ISBN 0-201-50255-0。

返回页首

附录 B：基本 HTM 例程

本节将介绍 HTM 例程。附带文档 [Szalay] 中含有用于每个例程的手册页，并且这些例程本身带有批注以支持 IntelliSense。

在下面，lat 和 lon 以十进制度数表示（南部和西部纬度为负），距离以海里（弧度分）表示。

HTM 库版本：fHtmVersion() 将返回版本字符串

例程将返回一个 nvarchar(max) 字符串来给出 HTM 库版本。

使用示例：

   print dbo.fHtmVersion()

返回的内容类似以下结果：

   'C# HTM.DLL V.1.0.0 1 August 2005'

生成 HTM 键：fHtmLatLon (lat, lon) 将返回 HtmID

例程将返回该 LatLon 点的 21 级深度的 HTM ID。

使用示例：

   update Place set HtmID = dbo.fHtmLatLon(Lat,Lon)

还有 fHtmXyz() 和 fHtmEq() 函数可供天文学家使用。

LatLon 到 XYZ：fHtmLatLonToXyz (lat,lon) 将返回点 (x, y, z)

例程将返回该 Lat Lon 点的笛卡尔坐标。

使用示例（这是标识函数）：

Select LatLon.lat, LatLon.lon-360
from fHtmLatLonToXyz(37.4,-122.4) as XYZ cross apply
fHtmXyzToLatLon(XYZ.x, XYZ.y, XYZ.z) as LatLon

还有 fHtmEqToXyz() 函数可供天文学家使用。

XYZ 到 LatLon：fHtmXyzToLatLon (x,y,z) 将返回点 (lat, lon)

例程将返回该 Lat Lon 点的笛卡尔坐标。

使用示例（这是标识函数）：

Select LatLon.lat, LatLon.lon-360
from fHtmLatLonToXyz(37.4,-122.4) as XYZ cross apply
fHtmXyzToLatLon(XYZ.x, XYZ.y, XYZ.z) as LatLon

还有 fHtmXyzToEq() 函数可供天文学家使用。

查看 HTM 键：fHtmToString (HtmID) 将返回 HTM 字符串

如果给定了 HtmID，则例程将返回一个 nvarchar(32)，其形式为 [N|S]t1t2t3...tn，其中每个三角数字 ti 均为 {0,1,2,3} 格式，用以说明三角网格中该深度的 HTM trixel。

使用示例：

print 'SQL Server development is at: ' +
dbo.fHtmToString(dbo.fHtmLatLon(47.646,-122.123))

返回结果：'N132130231002222332302'。

还有 fHtmXyz() 和 fHtmEq() 函数可供天文学家使用。

HTM trixel 中心点：fHtmToCenterpoint(HtmId) 将返回点 (x, y, z)

返回由 HtmID 指定的 HTM trixel 的笛卡尔中心点。

使用示例：

select * from fHtmToCenterPoint(dbo.fHtmLatLon(47.646,-122.123))

HTM trixel 角点：fHtmToCornerpoints(HtmId) 将返回点 (x, y, z)

返回由 HtmID 指定的 HTM trixel 的三个笛卡尔角点。

使用示例：

select * from fHtmToCornerPoints(dbo.fHtmLatLon(47.646,-122.123))

计算距离：fDistanceLatLon(lat1, lon1, lat2, lon2) 将返回距离

以海里（弧度分）为单位计算两点之间的距离。

使用示例：

declare @lat float, @lon float
select @lat = lat, @lon = lon
from  Place
where PlaceName = 'Baltimore' and State = 'MD'
select PlaceName,
dbo.fDistanceLatLon(@lat,@lon, lat,  lon) as distance
from Place

还有 fDistanceXyz() 和 fDistanceEq() 函数可供天文学家使用。

下面的例程可返回一个用作空间索引的表。所返回空间索引表的数据定义为：

SpatialIndexTable table (
HtmID   bigint   not null , -- HTM spatial key (based on lat/lon)
Lat     float    not null , -- Latitude in Decimal
Lon     float    not null , -- Longitude in Decimal
x       float    not null , -- Cartesian coordinates,
y       float    not null , -- derived from lat-lon
z       float    not null , --,
Type    char(1)  not null , -- place (P) or gauge (G)
ObjID   bigint   not null , -- object ID in table
distance float   not null , -- distance in arc minutes to object
primary key (HtmID, ObjID) )

查找邻近的对象：fHtmNearbyLatLon(type, lat, lon, radius) 将返回空间索引表

返回半径范围内特定类型的对象列表及它们到给定点的距离。该列表将按对象的位置由近到远排列。

使用示例：

   select distance, Place.*
from fHtmNearbyLatLon('P', 39.3, -76.6, 10) I join Place
on I.objID = Place.HtmID
order by distance

还有 fHtmGetNearbyEq() 和 fHtmGetNearbyXYZ() 函数可供天文学家使用。

查找最近的对象：fHtmNearestLatLon(type, lat, lon) 将返回空间索引表

返回包含距离该点最近的特定类型对象的列表。

使用示例：

   select distance, Place.*
from fHtmNearestLatLon('P', 39.3, -76.6) I join Place
on I.objID = Place.HtmID

还有 fHtmGetNearbyEq() 和 fHtmGetNearbyXYZ() 函数可供天文学家使用。

下面的例程将返回一个表，用以说明覆盖所需区域的一组 trixel（HIM 三角形）的 HtmIdStart 和 HtmIdEnd。表的定义为：

TrixelTable table (
HtmIdStart   bigint not null , -- min HtmID in trixel
HtmIdEnd     bigint not null   -- max HtmID in trixel
)

圆圈区域 HTM 覆盖：fHtmCoverCircleLatLon(lat, lon, radius) 将返回 trixel 表

返回覆盖指定圆圈的 trixel 表。

使用示例：

declare @answer nvarchar(max)
declare @lat float, @lon float
select @lat = lat, @lon = lon
from Place
where Place.PlaceName = 'Baltimore'
and State = 'MD'
set @answer = ' using fHtmCoverCircleLatLon() it finds:      '
select @answer = @answer
+ cast(P.placeName as varchar(max)) + ', '
+ str( dbo.fDistanceLatLon(@lat,@lon, I.lat, I.lon) ,4,2)
+ '  arc minutes distant.'
from SpatialIndex I join fHtmCoverCircleLatLon(@lat, @lon, 5)
On HtmID between HtmIdStart and HtmIdEnd -- coarse test
and type = 'P'                           -- it is a place
and dbo.fDistanceLatLon(@lat,@lon, lat, lon)
between 0.1 and 5  -- careful test
join Place P on I.objID = P.HtmID
order by dbo.fDistanceLatLon(@lat,@lon, I.lat, I.lon) asc
print 'The city within 5 arc minutes of Baltimore is: '
+ 'Lansdowne-Baltimore Highlands, 4.37 arc minutes away'

还有 fHtmCoverCircleEq() 函数可供天文学家使用。

HTM 覆盖的常规区域指定：fHtmCoverRegion(region) 将返回 trixel 表

返回覆盖指定区域的 trixel 表（本主题前面对区域进行了介绍）。

select S.*
from  (select ObjID
from fHtmCoverRegion('RECT LATLON 37 -109.55  41 -102.05')
loop join SpatialIndex
on HtmID between HtmIdStart and HtmIdEnd
and lat between 37 and 41
and lon between -109.05 and -102.048
and type = 'S') as G
join Station S on G.objID = S.StationNumber
OPTION (FORCE ORDER)

常规区域指定：fHtmRegionToNormalFormString(region) 将返回区域字符串

返回格式为 REGION {CONVEX {x y z d}* }* 的字符串，其中已从每个凸形删除多余的半空间；如 [Fekete] 中所述，凸形已被简化。

print dbo.fHtmToNormalForm('RECT LATLON 37 -109.55  41 -102.05')

下面的例程将返回一个表，用以说明覆盖所需区域的一组 trixel（HIM 三角形）的 HtmIdStart 和 HtmIdEnd。表的定义为：

RegionTable (convexID     bigint not null , -- ID of the convex, 0,1,...
halfSpaceID bigint not null     -- ID of the halfspace
-- within convex, 0,1,2,
x           float  not null     -- Cartesian coordinates of
y           float  not null     -- unit-normal-vector of
z           float  not null     -- halfspace plane
d           float  not null     -- displacement of halfspace
)                               -- along unit vector [-1..1]

将区域字符串转换为表：fHtmRegionToTable(region) 将返回区域表

返回一个表，用以将区域作为凸形组合来说明，其中每个凸形均为 x，y，z，d 半空间的交集。如 [Fekete] 中所述，凸形已被简化。本文第 4 节介绍了此函数的用法。

select *
from dbo.fHtmToNormalForm('RECT LATLON 37 -109.55  41 -102.05')

查找区域内的点：fHtmRegionObjects(region, type) 将返回对象表

返回一个表，其中包含空间索引中具有指定类型且位于区域内的对象的对象 ID。

select *            -- find Colorado places.
from Places join
where HtmID in
select objID
from dbo. fHtmRegionObjects('RECT LATLON 37 -109.55  41 -102.05','P')

常规区域诊断：fHtmRegionError(region ) 将返回消息

如果区域定义有效，则返回“OK”；否则，返回描述区域定义问题的诊断信息，并且后跟区域的语法定义。

print dbo.fHtmRegionError ('RECT LATLON 37 -109.55 41 -102.05')

转到原英文页面

返回页首

转载说明
作者：来自微软MSDN相关文章
网址：本文引用自http://msdn2.microsoft.com/zh-cn/library/aa964138.aspx