65.9K
CodeProject 正在变化。 阅读更多。
Home

使用 SQL Server 进行地理编码

作者：gstolarov
starIconstarIconstarIconstarIcon
emptyStarIcon
starIcon

4.68/5 (9投票s)

2010年3月1日

CPOL

11分钟阅读

viewsIcon

150278

downloadIcon

1620

使用 SQL Server 进行地理编码

引言

本文探讨了与地理编码相关的几个问题，特别是

  • 将 ESRI Shape 文件转换为 XML
  • 解析地址
  • 地理编码地址

所有的地理编码都实现为 SQL Server 函数，因此可以提供使用其他方法（如 Web 服务或带有地图软件的 COM 接口（即 MapPoint））无法达到的性能提升。

背景

我第一次尝试地理编码大概是在 3 年前。最初使用一些更广泛的可用服务，如 Google 或 Yahoo Web 服务，确实可以完成任务。但后来，我的实现开始遇到了瓶颈。namely，我有两个限制需要处理：

  • 将地址通过公共网络发送在某些我正在做的项目中是不可接受的。
  • 我偶尔需要对成千上万条遗留数据进行地理编码。

在网上搜索，我找到了 John Sample 的精彩文章 http://www.johnsample.com/articles/GeocodeWithSqlServer.aspx。他的研究非常有帮助。显然，你所需的一切都在那里，还有大量的免费数据——美国人口普查局的 TIGER/Line 数据，覆盖了美国每个县的每个街区的纬度/经度。一旦有了简单的文本文件，你只需将其加载到数据库中，然后使用插值技术，就可以相对轻松地找到美国每个地址的纬度/经度。

这对我来说已经奏效了几年。然后我开始发现一些问题，我主要将其归咎于 2005-2007 年的住房繁荣 :-)。看，这段时间显然有很多新建房屋，我开始注意到很多新地址在我开始时获得的 TIGER/Line 文件中并不存在。没关系，我想，我只需要获取最新的 TIGER/Line 文件，然后就可以继续了。好吧，猜猜怎么着，美国人口普查局不提供文本格式的 TIGER/Line 文件。现在他们只提供 ESRI 二进制 Shape 文件——这给我的计划带来了麻烦——我真的需要一种方法来获取最新的文本格式地理数据，以便将其加载到我的 SQL Server 中。

本文的其余部分描述了该项目的三个部分

  • 将 ESRI Shape 文件转换为 XML
  • 将 XML 数据加载到 SQL Server
  • SQL Server 函数，用于解析美国地址并进行地理编码

TIGER 文件

美国人口普查局提供拓扑集成地理编码和参照系统——TIGER - 数据。这是每个地区、街区、地标——几乎任何值得关注的事物的县级地理信息。特别是，我对可以按县下载的“所有线路”文件感兴趣（前往 http://www2.census.gov/cgi-bin/shapefiles2009/national-files，导航到你的州，然后是县，下载“所有线路”文件）。在这个文件中，对于每个街区的每个城市街区，你可以找到地址范围、街道名称、邮政编码以及当然的起始和结束街道街区的纬度/经度等信息。最初，此信息以 ASCII 文件下载的形式提供，但自 2006 年以来，此信息仅以 ESRI Shape 文件格式提供。

ESRI Shape 文件 / XML 转换

ESRI 似乎已成为空间数据存储和呈现的实际标准。至少他们的文件和文件支持变得越来越广泛。

也许提供 Shape 文件的高度概览是有价值的。Shape 文件实际上是一组包含至少 4 个文件

  • SHP - 包含实际形状的 Shape 文件——一个数据结构列表，主要包含对象类型和每个对象的点列表
  • SHX - Shape 索引文件，用于指示文件中每个对象的位置
  • DBF - 一个老式的 DBase 数据库文件，其中包含每个 Shape/对象的附加属性。对于 TIGER/Line，这将包含街道名称、邮政编码……但它也可以包含其他通用信息。
  • PRJ - 文件信息，例如文件中使用的编码/投影类型。

ESRI 是事实上的领导者这一事实并没有让我的工作轻松多少——将 ESRI Shape 文件转换为文本格式仍然花费了很长时间。最终，我没有尝试寻找一个可以为我完成这项工作的工具，而是开始寻找可以用来完成这项工作的库。

虽然我真的很想找到 C# 代码来完成工作，但最终我选择了 Shapefile C 库 - http://shapelib.maptools.org/，这是最简单且没有依赖项的。他们还有很多关于如何使其工作的示例，这让我的工作相对轻松——打开 DBF 和 SHP 文件，然后为其中的每个对象转储相应的 XML 节点。

对象类型

就地理编码而言，我只对折线感兴趣。并且为了简单起见，我只对直线插值感兴趣，所以即使是折线，我也只对起始点和终点感兴趣。因此，如果你的县区有很多长街区和弯曲的街道，你可能需要考虑让这个过程更精确一些。因此，为了减小文件大小，当 আমি 将对象导出到 XML 文件时，每当遇到折线对象时，我只会导出起始点和终点位置。

对于我参与过的一些其他项目，我需要导出多边形。对于这些，我导出所有点——所以当你查看代码时，不要感到惊讶。

州平面坐标转换

SHP 文件转换的这一部分不用于地址地理编码。但是，对于我参与过的一些其他项目，我不得不处理使用州平面坐标系统编码的 Shape 文件——所谓的兰勃特等角圆锥投影。冒着被地理专家嘲笑的风险，我对这种坐标系统的简单描述是我所理解的——为了提高单独小区域的测绘精度，一种地理编码类型是选择区域的中心，然后提供坐标而不是纬度/经度，而是类似以英尺为单位的向北/向西的坐标——（就像海盗地图——从老树 20 步向北，3 步向东）。听起来很傻，但它似乎相当普遍。至少我不得不处理它，所以我不得不构建一个模块，可以透明地检查 Shape 文件的投影系统并将其即时转换为纬度/经度。

作为转换的基础，我选择了蒙大拿州立图书馆 Montana State Library 发布的代码，并使其更通用，特别是添加了从 Shape 文件的 PRJ 文件中获取投影信息代码。

XML 文件格式

我选择创建的 XML 文件远非通用或可扩展，可能无法完全满足你的所有需求，但它确实为我完成了工作。

<layer>
<shp
        property='value' property='value' property='value'
        xFrom='lon' yFrom='lat' xTo='lon' yTo='lat'
>
      lon lat lon lat ...
</shp>
</layer>

在此文件中，<shp/> 节点为每个 Shape 对象重复。DBF 文件中的每个列/属性都成为 XML 文件中的一个属性（属性=值对）。

如果 Shape 对象是折线（想象一下城市街区），那么会添加 xFrom, yFrom, xTo, yTo 属性来指定折线的起始点和终点。对于任何其他对象类型，我会像 Google KML 文件一样添加经度/纬度点。

Shp2Xml 实用程序将假定所有 Shape 文件组件都在同一个文件夹中，并将输出转储到控制台。因此，要生成 XML 文件，请从命令行运行以下命令

 c:\>Shp2Xml SHP_file_name > xml_file_name

文件生成后，你可以使用 grep/findstr 之类的工具来过滤掉你不感兴趣的属性，以减小文件大小，因为这些文件可能会变得非常大——我们说的是数百兆字节。

XML 文件加载

在下载文件中，SQL 文件夹下有一个 LoadTiger.sql 文件，其中包含本文其余部分所述的所有代码。

加载部分非常简单

DECLARE @hDoc int
DECLARE @xml xml
set @xml=(select convert(xml,BulkColumn, 2)
From openrowset(Bulk 'C:\temp\geotrans\edges.xml', single_blob) [rowsetresults])
exec sp_xml_preparedocument @hDoc OUTPUT, @xml
SELECT *
into geoStreet
FROM OPENXML(@hDoc, 'layer/shp', 1)
WITH (
	tlid		int '@TLID',
	prfx		varchar(2)   ,
	street		varchar(50)  '@FULLNAME',
	type		varchar(4)   ,
	LFROMADD	varchar(20) '@LFROMADD',
	LTOADD		varchar(20) '@LTOADD',
	RFROMADD	varchar(20) '@RFROMADD',
	RTOADD		varchar(20) '@RTOADD',
	fraddr		int  ,
	toaddr		int  ,
	zipl		varchar(5)   '@ZIPL',
	zipr		varchar(5)   '@ZIPR',
	frlon		decimal(9,5) '@xFrom',
	frlat		decimal(9,5) '@yFrom',
	tolon		decimal(9,5) '@xTo',
	tolat		decimal(9,5) '@yTo'
)
exec sp_xml_removedocument @hDoc

此段将 c:\temp\geotrans\edges.xml （根据需要更改）加载到 geoStreet 表中。正如你所见，我只加载了以下列/属性：FULLNAME, LFROMADD, LTOADD, RFROMADD, RTOADD, ZIPL, ZIPR, xFrom, yFrom, xTo, yTo。

其余的加载代码用于数据清理。

delete from geoStreet 
where street is null 
	or lfromadd is null or ltoadd is null
	or rfromadd is null or rtoadd is null

此段从数据集中删除非街区元素——高速公路、上下匝道、河流……。另外，正如你所见，TIGER/Line 数据为街道的左侧/右侧提供了单独的信息。有时左侧和右侧的街道地址编号不如我希望的那样——要么一侧是非数字的，要么有时左侧地址是递增的，而右侧地址是递减的——很奇怪，是吧！为了使地理编码过程更容易，我希望从地址和到地址是可预测的递增顺序。因此，下面的代码会清理这部分。

update geoStreet set lFromAdd=rFromAdd where isnumeric(lFromAdd)=0
update geoStreet set rFromAdd=lFromAdd where isnumeric(rFromAdd)=0
update geoStreet set lToAdd=rToAdd where isnumeric(lToAdd)=0
update geoStreet set rToAdd=lToAdd where isnumeric(rToAdd)=0
update g set g.prfx=a.dir, g.type=a.tp, g.street=a.street
from geoStreet g cross apply fnParseAddr(g.street, 0) a

declare @b varchar(20), @tmpi int, @tmpn numeric(9,5), @tmpt numeric(9,5)
update geoStreet set 
	@b=LFROMADD, LFROMADD=LTOADD, LTOADD=@b
where   convert(int, LFROMADD)>convert(int, LTOADD)
	and convert(int, RFROMADD)<convert(int, RTOADD)
update geoStreet set 
	toaddr=(case when convert(int, RTOADD) > convert(int, LTOADD) _
	then RTOADD else LTOADD end),
	fraddr=(case when convert(int, RFROMADD) < convert(int, LFROMADD) _
	then RFROMADD else LFROMADD end)
where toaddr is null and fraddr is null
update geoStreet set 
	@tmpi=frAddr, frAddr=toAddr, toAddr=@tmpi,
	@tmpn=frlon, frlon=tolon, tolon=@tmpn,
	@tmpt=frlat, frlat=tolat, tolat=@tmpt
where frAddr > toAddr

唯一剩下的部分是处理 TIGER/Line 数据中的“所有线路”部分将街道名称作为单个字符串的事实。我决定手动将街道名称解析为方向/名称/类型元素，而不是尝试连接多个文件（稍后会详细介绍解析）。

update g set g.prfx=a.dir, g.type=a.tp, g.street=a.street
from geoStreet g cross apply fnParseAddr(g.street, 0) a

地址解析/地理编码

我发现需要处理的地址格式非常多，让我不知所措。有时我会将所有内容解析为单独的字段——号码、街道公寓等。有时我会获得地址行 1、城市和邮政编码。有时是地址行 1 和地址行 2。此外，有时我需要使用 C# 代码进行地理编码，其中存储过程易于执行，但有时我需要在 SQL Server 脚本中进行，其中存储过程也意味着使用游标以及循环、检查错误代码……的所有复杂性。所以我决定让自己轻松一些：我的地理编码将作为表值函数提供，该函数还将执行地址解析，并且只接受以下格式的单个地址字符串：

StrNum StrDir StrName StrType # apt, City State Zip

其中

  • StrNum 是数字
  • StrDir - N、S、E、W、North、South、East、West 之一
  • StrType - 可能的街道类型之一
  • # - 分隔可能的公寓号码
  • , - 第一个逗号将第一个地址行与第二个地址行分开
  • Zip - 5 位数字可选邮政编码
  • State - 2 位字符可选州
  • City - 第二个地址行上的其余所有内容

因此，fnParseAddr 的第一部分执行的正是这个——一次剥离字符串标记，试图解析地址。

查找第一个逗号以将地址拆分为第 1 行和第 2 行。

set @i = charindex(',', @str)
if (@i > 0) 
	select  @city = replace(ltrim(substring(@str, @i+1, 999)), ',', ' '),
		@str = replace(rtrim(substring(@str, 1, @i-1)), ',', ' ')

检查第一个标记是否为数字，如果是，则加载街道号码。

set @i = charindex(' ', @str)
if (@i > 0 and isnumeric(substring(@str, 1, @i))=1) 
	select  @num = rtrim(substring(@str, 1, @i-1)), 
		@str = ltrim(substring(@str, @i+1, 999))

检查第二个标记是否为街道方向（North/East/West/South）。

set @i = charindex(' ', @str)
if (@i > 0 and substring(@str, 1, @i-1) IN _
('N','No','North','South','S','West','W','East','E'))
	select  @dir = rtrim(substring(@str, 1, 1)), 
		@str = ltrim(substring(@str, @i+1, 999))

依此类推。

地址解析后，我们可以尝试确定它属于哪个街区。将所有街区加载到数据库中应该能使这项工作变得微不足道——只需一个简单的 select ，其中街道名称等于指定的内容，并且号码在范围内。但是，街道地址很少是正确的。大多数时候，人们不知道它是一条街还是一个大道。南北东西街道方向很少指定。邮政编码经常是错误的。因此，除了简单的匹配之外，我们还需要一个权重系统——每次匹配都会增加找到正确地址的可能性——如果邮政编码匹配，则得分加 3；如果地址在范围内，则加 2；如果街道类型或方向匹配，则加 1。使用最佳匹配记录。

这转化为以下 SQL 代码：

@wt = (case when @str=street then 5 else 0 end)
	+	(case when @zip IN (zipl, zipr) then 3 else 0 end)
	+ 	(case when @num between fraddr and toaddr then 2 else 0 end)
	+ 	(case when @dir=prfx then 1 else 0 end)
	+ 	(case when @tp=type then 1 else 0 end)

现在已经找到了最匹配我们地址的记录，我们需要找到实际的纬度/经度。为此，我们使用插值：如果一个街区包含地址从 100 到 200，而我们正在寻找房子号 125，那么它应该在街区开始处的大约四分之一处。

select @ratio = 0.5
if @num between @fraddr and @toaddr or @num between @toaddr and @fraddr 
	select @ratio = 1.0 *(@num - @fraddr) / (@toaddr - @fraddr)
select 	@lat = (@ratio * (@tolat - @frlat) + @frlat) , 
	@lon = (@ratio * (@tolon - @frlon) + @frlon)

Using the Code

运行 LoadTiger.sql 脚本后，它将创建 geoStreet 表和 fnParseAddr 函数。该表由函数用于地理编码地址。该函数接受两个参数：

  • Address - 有关接受的地址格式，请参阅上一节。
  • ifGeocode - 0/1 标志，用于指定是只需要地址解析还是需要执行完整的地理编码。

该函数返回一个仅包含 1 行的表：

num varchar(20) 街道号码
dir varchar(1) 街道方向
street varchar(100) 街道名称
tp varchar(20) 街道类型
apt varchar(20) 公寓
city varchar(100) city
st varchar(2) 状态
zip varchar(5) zip
lat decimal(9,5) 纬度
lon decimal(9,5) 经度

要仅解析地址，请使用以下 SQL 命令：

select * from fnParseAddr('123 main st', 0)

要解析和地理编码地址，请发出以下 select 语句：

select * from fnParseAddr('123 main st', 1)

要对数据库中的某个表执行批量地理编码，请发出以下命令：

select tbl.stname, a.*
from SomeTable tbl 
	cross apply fnParseAddr(tbl.stname, 1) a

或者，类似地更新信息：

update tbl set tbl.lat=a.lat, tbl.lon=a.lon 
from SomeTable tbl 
	cross apply fnParseAddr(tbl.stname, 1) a
where a.lat != 0

注释

项目的某些部分可能受益于改进。以下是我希望在某个时候解决的问题列表：

  • 更灵活的地址解析。目前，对于像 123 Avenue 52 E 这样的地址，我会遇到问题——它不太符合常规地址格式。另外，解析公寓号码（目前用 # 分隔）还有改进的空间。
  • 支持除纬度/经度或州平面以外的坐标投影系统。
  • 性能改进——虽然我在双核 3GHz Xeon 和 3GB RAM 上获得了令人满意的 1,700 条记录/秒的性能，但通过创建托管代码来解析地址，性能可能会进一步提高。
© . All rights reserved.