Using Table Valued Functions in SQL Server 2005 to Implement a Spatial Data Library
Gyorgy Fekete and Alex Szalay
Johns Hopkins University
Jim Gray
Microsoft (contact author)
November 2005
Applies to:
Microsoft SQL Server 2005
Summary: This article explains how to add spatial search functions (point-near-point and point in polygon) to Microsoft SQL Server 2005 using C# and table-valued functions. It is possible to use this library to add spatial search to your application without writing any special code. The library implements the public-domain C# Hierarchical Triangular Mesh (HTM) algorithms from Johns Hopkins University. That C# library is connected to SQL Server 2005 using a set of scalar-valued and table-valued functions. These functions act as a spatial index. (26 printed pages)
The sample code for this article is included with SQL Server 2005. It is the spatial CLR sample, and by default it installs to drive:\Program Files\Microsoft SQL Server\90\Samples\Engine\Programmability\CLR\Spatial.
Contents
Introduction
Table Valued Functions: The Key Idea
Using Table-Valued Functions to Add a Spatial Index
The Datasets
Typical Queries
Conclusion
Appendix A: References
Appendix B: Basic HTM Routines
Introduction
Spatial data searches are common in both commercial and scientific applications. We developed a spatial search system in conjunction with our effort to build the SkyServer (http://skyserver.sdss.org) for the astronomy community. The SkyServer is a multi-terabyte database that catalogs about 300 million celestial objects. Many of the questions astronomers want to ask of it involve spatial searches. Typical queries include, "What is near this point," "What objects are inside this area," and "What areas overlap this area?"
For this article, we have added the latitude/longitude (lat/lon) terrestrial sphere (the earth) grid to the astronomer's right ascension/declination (ra/dec) celestial sphere (the sky) grid. The two grids have a lot in common, but the correspondence is not exact; the traditional order lat-lon corresponds to dec-ra. This order reversal forces us to be explicit about the coordinate system. We call the Greenwich-Meridian-Equatorial terrestrial coordinate system the LatLon coordinate system. The library supports three coordinate systems:
- Greenwich Latitude-Longitude, called LatLon.
- Astronomical right-ascension—declination called J2000.
- Cartesian (x, y, z) called Cartesian.
Astronomers use arc minutes as their standard distance metric. A nautical mile is an arc minute, so the distance translation is very natural. Many other concepts are quite similar. To demonstrate these, this article will show you how to use the spatial library to build a spatial index on two USGS datasets: US cities and US stream-flow gauges. Using these indexes and some spatial functions, the article provides examples of how to search for cities near a point, how to find stream gauges near a city, and how to find stream gauges or cities within a state (polygonal area).
We believe this approach is generic. The spatial data spine schema and spatial data functions can be added to almost any application to allow spatial queries. The ideas also apply to other multi-dimensional indexing schemes. For example, the techniques would work for searching color space or any other low-dimension metric space.
Table Valued Functions: The Key Idea
The key concept of relational algebra is that every relational operator consumes one or more relations and produces an output relation. SQL is syntactic sugar for this idea, allowing you to define relations (data definition language) and to manipulate relations with a select-insert-update-delete syntax.
Defining your own scalar functions lets you make some extensions to the relational database—you can send mail messages, you can execute command scripts, and you can compute non-standard scalars and aggregate values such as tax() or median().
However, if you can create tables, then you can become part of the relational engine—both a producer and consumer of relational tables. This was the idea of OLEDB, which allows any data source to produce a data stream. It is also the idea behind the SQL Server 2000 Table Valued Functions.
Implementing table valued functions in Transact-SQL is easy:
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
This is fine if your function can be done entirely in Transact-SQL. But implementing OLEDB data sources or Table Valued Functions outside of Transact-SQL is a real challenge in SQL Server 2000.
The common language runtime (CLR) integration of SQL Server 2005 makes it easy to create a table-valued function. You create a list, an array, or any IEnumerable object (anything you can do foreach on), and then you cast it as a table.
[SqlFunction(TableDefinition = "x float, y float" ,
FillRowMethodName = "FillPair")]
public static IEnumerable csPoints( )
{
int[,] points = { { 1, 2 }, { 3, 4 } };
return (IEnumerable) points;
}
You compile this in Visual Studio and click Deploy. The table-valued function is installed in the database.
Using Table-Valued Functions to Add a Spatial Index
There is a lot of confusion about indexes. Indexes are really simple—they are tables with a few special properties:
- SQL Server has only one kind of associative (by value) index—a B-tree. The B-tree can have multi-field keys, but the first field carries most of the selectivity.
- Conceptually, the B-Tree index is a table consisting of the B-Tree key fields, the base table key fields, and any included fields that you want to add to the index.
- B-tree indexes are sorted according to the index key, such as ZIP code or customer ID, so that lookup or sequential scan by that key is fast.
- Indexes are often smaller than the base table, carrying only the most important attributes, so that looking in the index involves many fewer bytes than examining the whole table. Often, the index is so much smaller that it can fit in main memory, thereby saving even more disk accesses.
- When you think you are doing an index lookup, you are either searching the index alone (a vertical partition of the base table), or you are searching the index, and then joining the qualifying index rows to rows in the base table using the base-table primary key (a bookmark lookup).
The central idea is that the spatial index gives you a small subset of the data. The index tells you where to look and often carries some helpful search information with it (called included columns or covering columns by the experts.) The selectivity of an index tells how big this initial reduction is (the coarse subset of Figure 1). When the subset is located, a careful test examines each member of the subset and discards false positives. That process is indicated by the diamond in Figure 1. A good index has few false positives. The metaphor shown in Figure 1 (the coarse subset and the careful test) is used throughout this article.
.gif "Aa964138.tvf_figure_1(en-US,SQL.90).gif")
Figure 1
B-trees and table-valued functions can be combined as follows to let you build your own spatial index that produces coarse subsets:
- Create a function that generates keys that cluster related data together. For example, if items A and B are related, then the keys for A and B should be nearby in the B-tree key space.
- Create a table-valued function that, given a description of the subset of interest, returns a list of key ranges (a "cover") containing all the pertinent values.
You cannot always get every key to be near all its relatives because keys are sorted in one dimension and relatives are near in two-dimensional space or higher. However, you can come close. The ratio of false-positives to correct answers is a measure of how well you are doing.
The standard approach is to find some space filling curve and thread the key space along that curve. Using the standard Mercator map, for example, you can assign everyone in the Northwest to the Northwest key range, and assign everyone in the Southeast to the Southeast key range. Figure 2 shows the second order space-filling curve that traverses all these quadrants, assigning keys in sequence. Everyone in the Northwest-Southwest quadrant has the key prefix nwsw. If you have an area like the circle shown in Figure 2, you can look in the key range
key between 'nwsw' and 'nwse'
.gif "Aa964138.tvf_figure_2(en-US,SQL.90).gif")
Figure 2
This search space is eight times smaller than the whole table and has about 75 percent false positives (indicated by the area outside the circle but inside the two boxes). This is not a great improvement, but it conveys the idea. A better index would use a finer cell division. With fine enough cells, the converging area could have very few false positives. A detailed review of space-filling curves and space-partitioning trees can be found in the books of Hanan Samet [Samet].
Now we are going to define a space-filling curve—the Hierarchical Triangular Mesh (HTM) that works particularly well on the sphere. The earth is round and the celestial sphere is round, so this spherical system is very convenient for geographers and astronomers. We could do similar things for any metric space. The space-filling curve gives keys that are the basis of the spatial index. Then, when someone has a region of interest, our table valued function will give them a good set of key-ranges to look at (the coarse filter of Figure 1). These key ranges will cover the region with spherical triangles, called trixels, much as the two boxes in Figure 2 cover the circle. The search function need only look at all the objects in the key ranges of these trixels to see if they qualify (the careful test in Figure 1).
To make this concrete, assume we have a table of Objects
create table Object (objID bigint primary key,
lat float, -- Latitude
lon float, -- Longitude
HtmID bigint) -- The HTM key
and a distance function dbo.fDistanceLatLon(lat1, lon1, lat2, lon2) that gives the distance in nautical miles (arc minutes) between two points. Further assume that the following table-valued function gives us the list of key ranges for HtmID points that are within a certain radius of a lat-lon point.
define function
fHtmCoverCircleLatLon(@lat float, @lon float, @radius float)
returns @TrixelTable table(HtmIdStart bigint, HtmIdEnd bigint)
Then the following query finds points within 40 nautical miles of San Francisco (lat,lon) = (37.8,-122.4):
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
We now must define the HTM key generation function, the distance function, and the HTM cover function. That's what we do next using two United States Geological spatial datasets as an example. If you are skeptical that this scales to billions of objects, go to http://skyserver.sdss.org and look around the site. That Web site uses this same code to do its spatial lookup on a multi-terabyte astronomy database.
This article is about how you use SQL Table Valued Functions and a space-filling curve like the HTM to build a spatial index. As such, we treat the HTM code itself as a black box documented elsewhere [Szalay], and we focus on how to adapt it to our needs within an SQL application.
The Datasets
The US Geological Survey gathers and publishes data about the United States. Figure 3 shows the locations of 18,000 USGS-maintained stream gauges that measure river water flows and levels. The USGS also publishes a list of 23,000 place names and their populations.
.gif "Aa964138.tvf_figure_3(en-US,SQL.90).gif")
Figure 3
USGS Populated Places (23,000 Cities)
The USGS published a list of place names and some of their attributes in 1993. There are newer lists at the USGS Web site but they are fragmented by state, so it is difficult to get a nationwide list. The old list will suffice to demonstrate spatial indicies. The data has the following format:
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
)
To speed name lookups, we add a name index, but the data is clustered by the spatial key. Nearby objects are co-located in the clustering B-tree and thus on the same or nearby disk pages.
create index Place_Name on Place(PlaceName)
All except the HtmID data can be downloaded from the USGS Web site. The SQL Server 2005 data import wizard can be used to import the data (we have already done that in the sample database). The HtmID field is computed from the Lat Lon by:
update Place set HtmID = dbo.fHtmLatLon(lat, lon)
USGS Stream Gauges (17,000 Instruments)
The USGS has been maintaining records of river flows since 1854. As of 1 Jan 2000, they had accumulated over 430 thousand years of measurement data. About six thousand active stations were active, and about four thousand were online. The gauges are described in detail at http://waterdata.usgs.gov/nwis/rt. An NOAA site shows the data from a few hundred of the most popular stations in a very convenient way:http://weather.gov/rivers_tab.php.
Our database has just the stations in the continental United States (see Figure 3). There are also stations in Guam, Alaska, Hawaii, Puerto Rico, and the Virgin Islands that are not included in this database. The stream gauge station table is:
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) )
As before, the HtmID field is computed from the Lat Lon fields by:
update Station set HtmID = dbo.fHtmLatLon(lat, lon)
There are up to 18 stations at one location, so the primary key must include the station number to make it unique. However, the HTM key clusters all the nearby stations together in the B-tree. To speed lookups, we add a station number and a name index:
create index Station_Name on Station(StationName)
create index Station_Number on Station(StationNumber)
The Spatial Index Table
Now we are ready to create our spatial index. We could have added the fields to the base tables, but to make the stored procedures work for many different tables, we found it convenient to just mix all the objects together in one spatial index. You could choose (type,HtmID) as the key to segregate the different types of objects; we chose (HtmID, key) as the key so that nearby objects of all types (cities and steam gagues) are clustered together. The spatial index is:
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) )
The Cartesian coordinates will be explained later in this topic. For now, it is enough to say that the function fHtmCenterPoint(HtmID) returns the Cartesian (x,y,z) unit vector for the centerpoint of that HTM triangle. This is the limit point of the HTM, as the center is subdivided to infinitely small trixels.
The SpatialIndex table is populated from the Place and Station tables as follows:
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
To clean up the database, we execute:
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
A Digression: Cartesian Coordinates
You can skip this if you like. It is not needed to use the library. The HTM code heavily uses a trick to avoid spherical geometry: it moves from the 2D surface of the sphere to 3D. This allows very quick tests for "inside a polygon" and for "nearby a point" queries.
Every lat/lon point on the sphere can be represented by a unit vector in three-dimensional space v = (x,y,z). The north and south poles (90° and -90°) are v = (0,0,1), and v = (0,0,-1) respectively. Z represents the axis of rotation, and the XZ plane represenst the Prime (Greenwich) Meridian, having longitude 0° or longitude 180°. The formal definitions are:
x = cos(lat)cos(lon)
y =cos(lat)sin(lon)
z = sin(lat)
These Cartesian coordiates are used as follows: given two points on the unit sphere, p1=(x1,y1,z1) and p2 = (x2,y2,z2), then their dot product, p1•p2 = x1*x2+y1*y2+z1*z2, is the cosine of the angle between these two points. It is a distance metric. Figure 4 shows how Cartesian coordinates allow quick tests for point-in-polygon and point-near-point. Each lat/lon point has a corresponding (x,y,z) unit vector.
.gif "Aa964138.tvf_figure_4(en-US,SQL.90).gif")
Figure 4
If we are looking for points within 45 nautical miles (arc minutes) of point p1, that is at most 45/60 degrees away from p1. The dot product of such points with p1 will be less than d=acos(radians(45/60). The "is nearby" test becomes { p2 | p2•p1 < d}, which is a very quick test. In Figure 5, each great or small circle is the intersection of a plane with the circle. A point is inside the circle if its dot product with the plane's normal vector is less than cos(è) where 2è is the circle's arc-angle diameter.
.gif "Aa964138.tvf_figure_5(en-US,SQL.90).gif")
Figure 5
Cartesian coordinates also allow a quick test for point-inside-polygon. All our polygons have great-circle or small-circle edges. Such edges lie along a plane intersecting the sphere. Therefore, the edges can be defined by the unit vector, v, normal to the plane and by a shift along that vector. For example, the equator is the vector v = (0,0,1) and shift zero. Latitude 60° is defined by vector v = (0,0,1) with a shift of 0.5, and a 60° circle around Baltimore is defined by vector v = (0.179195, -0.752798, 0.633392) with a shift of 0.5. A place, p2, is within 60° of Baltimore if p2•v < 0.5. The same idea lets us decide if a point is inside or outside a HTM triangle by evaluating three such dot products. That is one of the main reasons the HTM code is so efficient and fast.
We have implemented several helper procedures to convert from LatLon to Cartesian coordiantes:
fHtmXyz(HtmID) returns the xyz vector of the centerpoint of an HtmID.
fHtmLatLonToXyz(lat,lon) returns an xyz vector.
fHtmXyzToLatLon(x,y,z) returns a lat,lon vector.
They are used later in the article and documented in the the API spec and Intellisense [Fekete].
The library here defaults to 21-deep HTM keys (the first level divides the sphere into eight faces and each subsequent level divides the speherical triangle into four sub-triangles.) Table 1 indicates that a 21-deep trixel is fairly small. The code can be modified to go 31-deep before the 64-bit representation runs out of bits.
In Table 1, each HTM level subdivdes the sphere. For each level, this table shows the area in square degrees, arc minutes, arc seconds, and meters. The Trixel column shows some charactic sizes: the default 21-deep trixels is about .3 arc second2. The USGS data has about half an object per 12-deep trixel.
Table 1
| HTM depth | Area | objects / trixel | |||||
|---|---|---|---|---|---|---|---|
| deg2 | arc min2 | arc sec2 | earth m2 | trixel | SDSS | USGS | |
| sphere | 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 | |||
Typical Queries
You should now be ready to do a few queries.
Query 1: Find Points Near Point: Find Towns Near A Place
The most common query is to find all places nearby a certain place or point. Consider the query, "Find all towns within 100 nautical miles of Baltimore, Maryland." The HTM triangles covering a 100 nautical mile circle (100 arc minutes from) Baltimore are obtained by
select * -- find a HTM cover 100 NM around Baltimore
from fHtmCoverCircleLatLon(39.3, -76.6, 100)
This returns the Trixel Table shown in Table 2. That is, the fHtmCoverCircleLatLon() function returns a set of HTM triangles that "cover" the circle (in this case, a single trixel). The HTM keys of all objects inside the circle are also inside one of these triangles. Now we need to look in all those triangles and discard the false positives (the careful test of Figure 1). We will order the answer set by the distance from Baltimore, so that if we want the closest place, we can just select the TOP 1 WHERE distance > 0 (we want to exclude Baltimore itself from being closest).
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
Table 2. The Baltimore circle HTM cover
| HtmIdStart | HtmIdEnd |
|---|---|
| 14023068221440 | 14027363188735 |
The cover join returns 2,928 rows (the coarse test); 1,122 of them are within 100 air miles (the careful test). This gives us 61 percent false positives—all within 9 milliseconds.
These are such common tasks that there are standard functions for them:
fHtmNearbyLatLon(type, lat, lon, radius)
fHtmNearestLatLon(type, lat, lon)
so the preceding query becomes
select ObjID, distance
from fHtmNearestLatLon('P', 39.3, -76.61)
Query 2: Find Places Inside a Box
Applications often want to find all the objects inside a square view-port when displaying a square map or window. Colorado is almost exactly square with corner points (41°N-109°3'W) in the NW corner and (37°N-102° 3'E) in the SW corner. The state's center point is (39°N, -105°33'E) so one can cover that square with a circle centered at that point.
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)
This example returns 1,030 stream gauges in about 46 milliseconds. Five other Colorado gauges are right on the border that wanders south of 37° by up to one nautical mile. These extra five stations appear when the southern latitude is adjusted from 37° to 36.98°. (GIS systems and astronomical applications often want a buffer zone around a region. The HTM code includes support for buffer zones, and they are much used in real applications. Look at reference [Szalay] to see how this is done.) The cover circle returns 36 triangles. The join with the SpatialIndex table returns 1,975 gauges. That's 47 percent false positives. The next section shows how to improve on this by using the HTM regions to specify a polygon cover rather than a cover for a circle.
The FORCE ORDER clause is an embarrassment—if missing, the query runs ten times longer because the optimizer does a nested-loops join with the spatial index as the outer table. Perhaps if the tables were larger (millions of rows), the optimizer would pick a different plan, but we cannot count on that. Paradoxically, the optimizer chose the correct plan without any hints for all the queries in the previous section.
Query 3: Find Places Inside a Polygon
The HTM code lets us specify the area as a circle, a rectangle, a convex hull, or a union of these regions. In particular, the HTM library allows you to specify a region using the following linear syntax:
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
To give examples of region specifications:
CIRCLE. A point specification and a 1.75 nautical mile (arc minute) radius.
'CIRCLE LATLON 39.3 -76.61 100' 'CIRCLE CARTESIAN 0.1792 -0.7528 0.6334 100'RECT. Two corner points defining the minimum and maximum of the lat, lon. The longitude coordinates are interpreted in the wrap-around sense, that is, lonmin=358.0 and lonmax=2.0, is a four-degree-wide range. The latitudes must be between the North and South Poles. The rectangle edges are constant latitude and longitude lines, rather than the great-circle edges of CHULL and CONVEX.
'RECT LATLON 37 -109.55 41 -102.05'CHULL. Three or more point specifications define a spherical convex hull with edges of the convex hull connecting adjacent points by great circles. The points must be in a single hemisphere, otherwise an error is returned. The order of the points is irrelevant.
'CHULL LATLON 37 -109.55 41 -109.55 41 -102.051 37 -102.05'CONVEX. Any number (including zero) of constraints in the form of a Cartesian vector (x,y,z) and a fraction of the unit length of the vector.
'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. A region is the union of zero or more circle, rect, chull, and convex areas.
'REGION CONVEX 0.7 0.7 0.0 -0.5 CIRCLE LATLON 18.2 -22.4 1.75'
Any of these region descriptions can be fed to the fHtmCoverRegion() routine that returns a trixel table describing a set of trixels (triangular areas) covering that region. The simpler code for the Colorado query is:
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)
This unusual query format is required to tell the optimizer exactly the order in which to perform the join (to make the "force order" option work correctly). It is difficult to modify the optimizer in this way, but until table-valued functions have statistics, they are estimated to be very expensive. You have to force them into the inner loop join.
The query returns 1,030 stream gauges and has 1,365 candidates from the cover, so there are 25 percent false positives. Note that the rectangle cover is better than the circular cover, which had 61 percent false positives. There is polygon syntax for non-rectangular states, but this article is about table valued functions, not about the HTM algorithms. You can see the HTM code in the project, and also in the documentation for the project.
A similar query can be cast as a convex hull.
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)
The query returns 1,030 stream gauges and has 1,193 candidates from the cover, so there are 14 percent false positives. The convex hull cover is even better than the equivalent rectangular cover in this case.
Query 4: Advanced Topics—Complex Regions
The previous examples gave the syntax for regions and a discussion of point-near-point and point-in-rectangle searches. Regions can get quite complex. They are Boolean combinations of convex areas. We do not have the space here to explain regions in detail, but the HTM library in the accompanying project has the logic to do Boolean combinations of regions, simplify regions, compute region corner points, compute region areas, and also has many other features. Those ideas are described in [Fekete], [Gray], and [Szalay].
To give a hint of these ideas, consider the state of Utah. Its boundaries are approximately defined by the union of two rectangles:
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.
Now we can find all stream gauges in Utah with the query:
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)
The cover returns 38 trixels. The join returns 775 stations. The careful test finds 670 stations in Utah, and two Wyoming stations that are right on the border (14 percent false positives).
Most states require much more complex regions. For example, a region string to approximate California is:
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)
The cover returns 108 trixels, which cover 2,132 stations. Of these, 1,928 are inside California, so the false positives are about 5 percent—but the careful test is nontrivial.
That same query, done for places rather than stations, with the careful test included, looks like this:
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)
This uses the convex-halfspace representation of California and the techniques described in [Gray] to quickly test if a point is inside the California convex hull. It returns 885 places, 7 of which are on the Arizona border with California (the polygon approximates California). It runs in 0.249 seconds on a 1GHz processor. If you leave off the "OPTION(FORCE ORDER)" clause it runs slower, taking 247 seconds.
Because this is such a common requirement, and because the code is so tricky, we added the procedure fHtmRegionObjects(Region,Type) that returns object IDs from SpatialIndex. This procedure encapsulates the previously shown tricky code, so the two California queries become:
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'))
The Colorado and Utah queries are also simplified by using this routine.
Conclusion
The HTM spatial indexing library presented here is interesting and useful in its own right. It is a convenient way to index data for point-in-polygon queries on the sphere. But, the library is also a good example of how SQL Server and other database systems can be extended by adding a class library that does substantial computation in a language like C#, C++, Visual Basic, or Java. The ability to implement powerful table-valued functions and scalar functions and integrate these queries and their persistent data into the database is a very powerful extension mechanism that starts to deliver on the promise of Object-Relational databases. This is just a first step. In the next decade, programming languages and database query languages are likely to get even better data integration. This will be a boon to application developers.
For more information:
https://msdn.microsoft.com/sql/
Project Editor: Susanne Bonney
Appendix A: References
- [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, April 2004
- [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, May 15, 2005, 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.
Appendix B: Basic HTM Routines
This section describes the HTM routines. The companion document [Szalay] has a manual page for each routine, and the routines themselves are annotated to support IntelliSense.
In what follows, lat and lon are in decimal degrees (southern and western latitudes are negative), and distances are in nautical miles (arc minutes.)
HTM library version: fHtmVersion() returns versionString
The routine returns an nvarchar(max) string giving the HTM library version.
Example use:
print dbo.fHtmVersion()
Returns something like:
'C# HTM.DLL V.1.0.0 1 August 2005'
Generating HTM keys: fHtmLatLon (lat, lon) returns HtmID
The routine returns the 21-deep HTM ID of that LatLon point.
Example use:
update Place set HtmID = dbo.fHtmLatLon(Lat,Lon)
There are also fHtmXyz() and fHtmEq() functions for astronomers.
LatLon to XYZ: fHtmLatLonToXyz (lat,lon) returns Point (x, y, z)
The routine returns the Cartesian coordinates of that Lat Lon point.
Example use (this is the identity function):
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
There is also an fHtmEqToXyz() function for astronomers.
XYZ to LatLon: fHtmXyzToLatLon (x,y,z) returns Point (lat, lon)
The routine returns the Cartesian coordinates of that Lat Lon point.
Example use (this is the identity function):
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
There is also an fHtmXyzToEq() function for astronomers.
Viewing HTM keys: fHtmToString (HtmID) returns HtmString
Given an HtmID, the routine returns a nvarchar(32) in the form [N|S]t1t2t3...tn, where each triangle number ti is in {0,1,2,3} describing the HTM trixel at that depth of the triangular mesh.
Example use:
print 'SQL Server development is at: ' +
dbo.fHtmToString(dbo.fHtmLatLon(47.646,-122.123))
which returns: 'N132130231002222332302'.
There are also fHtmXyz() and fHtmEq() functions for astronomers.
HTM trixel Centerpoint: fHtmToCenterpoint(HtmId) returns Point (x, y, z)
Returns the Cartesian center point of the HTM trixel specified by the HtmID.
Example use:
select * from fHtmToCenterPoint(dbo.fHtmLatLon(47.646,-122.123))
HTM trixel corner points: fHtmToCornerpoints(HtmId) returns Point (x, y, z)
Returns the three Cartesian corner points of the HTM trixel specified by the HtmID.
Example use:
select * from fHtmToCornerPoints(dbo.fHtmLatLon(47.646,-122.123))
Computing distances: fDistanceLatLon(lat1, lon1, lat2, lon2) returns distance
Computes the distance, in nautical miles (arc minutes) between two points.
Example use:
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
There are also fDistanceXyz() and fDistanceEq() functions for astronomers.
The following routines return a table which serves as a spatial index. The returned spatial index table has the data definition:
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) )
Finding nearby objects: fHtmNearbyLatLon(type, lat, lon, radius) returns SpatialIndexTable
Returns a list of objects within the radius distance of the given type and their distance from the given point. The list is sorted by nearest object.
Example use:
select distance, Place.*
from fHtmNearbyLatLon('P', 39.3, -76.6, 10) I join Place
on I.objID = Place.HtmID
order by distance
There are also fHtmGetNearbyEq () and fHtmGetNearbyXYZ() functions for astronomers.
Finding the nearest object: fHtmNearestLatLon(type, lat, lon) returns SpatialIndexTable
Returns a list containing the nearest object of the given type to that point.
Example use:
select distance, Place.*
from fHtmNearestLatLon('P', 39.3, -76.6) I join Place
on I.objID = Place.HtmID
There are also fHtmGetNearestEq() and fHtmGetNearestXYZ() functions for astronomers.
The following routines return a table describing the HtmIdStart and HtmIdEnd of a set of trixels (HTM triangles) covering the area of interest. The table definition is:
TrixelTable table (
HtmIdStart bigint not null , -- min HtmID in trixel
HtmIdEnd bigint not null -- max HtmID in trixel
)
Circular region HTM cover: fHtmCoverCircleLatLon(lat, lon, radius) returns trixelTable
Returns a trixel table covering the designated circle.
Example use:
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'
There is also an fHtmCoverCircleEq() function for astronomers.
General region specification to HTM cover: fHtmCoverRegion(region) returns trixelTable
Returns a trixel table covering the designated region (regions are described earlier in this topic).
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)
General region simplification: fHtmRegionToNormalFormString(region) returns regionString
Returns a string of the form REGION {CONVEX {x y z d}* }* where redundant halfspaces have been removed from each convex; the convex has been simplified as described in [Fekete].
print dbo.fHtmToNormalForm('RECT LATLON 37 -109.55 41 -102.05')
The following routine returns a table describing the HtmIdStart and HtmIdEnd of a set of trixels (HTM triangles) covering the area of interest. The table definition is:
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]
Cast RegionString as Table: fHtmRegionToTable(region) returns RegionTable
Returns a table describing the region as a union of convexes, where each convex is the intersection of the x,y,z,d halfspaces. The convexes have been simplified as described in [Fekete]. Section 4 of this article describes the use of this function.
select *
from dbo.fHtmToNormalForm('RECT LATLON 37 -109.55 41 -102.05')
Find Points Inside a Region: fHtmRegionObjects(region, type) returns ObjectTable
Returns a table containing the object IDs of objects in SpatialIndex that have the designated type and are inside the region.
select * -- find Colorado places.
from Places join
where HtmID in
select objID
from dbo. fHtmRegionObjects('RECT LATLON 37 -109.55 41 -102.05','P')
General region diagnostic: fHtmRegionError(region ) returns message
Returns "OK" if region definition is valid; otherwise, returns a diagnostic describing what is wrong with the region definition, followed by a syntax definition of regions.
print dbo.fHtmRegionError ('RECT LATLON 37 -109.55 41 -102.05')