Geocoding For Data Scientists - KDnuggets

由柏拉图重新发布

关注： 0

当数据科学家需要了解有关数据“位置”的所有信息时，他们通常会求助于地理信息系统 (GIS)。 GIS 是一套复杂的技术和程序，具有多种用途，但华盛顿大学提供了相当全面的定义，表示“地理信息系统是关联或连接的事物或对象的复杂排列，其目的是交流有关地球表面特征的知识”（Lawler 等人）。 GIS 包含处理空间数据（从采集到可视化）的广泛技术，即使您不是 GIS 专家，其中许多技术也是有价值的工具。本文全面概述了地理编码，并用 Python 演示了几个实际应用程序。具体来说，您将使用地址确定纽约州纽约市一家披萨店的确切位置，并将其连接到附近公园的数据。虽然演示使用 Python 代码，但核心概念可以应用于许多编程环境，以将地理编码集成到您的工作流程中。这些工具为将数据转换为空间数据提供了基础，并为更复杂的地理分析打开了大门。

XXXXX

地理编码最常被定义为将地址数据转换为地图坐标。通常，这涉及检测地址中的街道名称，将该街道与数据库中现实世界对应街道的边界进行匹配，然后使用街道号码估计在街道上放置该地址的位置。作为示例，让我们对纽约百老汇的一家披萨店地址进行简单的手动地理编码：2709 Broadway, New York, NY 10025。第一个任务是为该位置的道路系统找到适当的 shapefile您的地址。请注意，在本例中，地址的城市和州是“纽约州纽约市”。幸运的是，纽约市公布了详细的道路信息纽约市开放数据页面（CSCL PUB）。其次，检查街道名称“百老汇”。您现在知道该地址可以位于纽约市任何名为“Broadway”的街道上，因此您可以执行以下 Python 代码来查询 NYC Open Data SODA API 以查找所有名为“Broadway”的街道。

import geopandas as gpd
import requests
from io import BytesIO # Request the data from the SODA API
req = requests.get( "https://data.cityofnewyork.us/resource/gdww-crzy.geojson?stname_lab=BROADWAY"
)
# Convert to a stream of bytes
reqstrm = BytesIO(req.content)
# Read the stream as a GeoDataFrame
ny_streets = gpd.read_file(reqstrm)

此查询有 700 多个结果，但这并不意味着您必须检查 700 条街道才能找到您的披萨。将数据可视化，您可以看到有 3 条主要的百老汇街道和一些较小的街道。

XXXXX

原因是每条街道都被分成大致对应于一个街区的部分，从而可以更精细地查看数据。该过程的下一步是使用邮政编码和街道号码准确确定地址位于这些部分的哪个部分。数据集中的每个街道段都包含街道左侧和右侧建筑物地址的地址范围。同样，每个路段都包含街道左侧和右侧的邮政编码。为了找到正确的路段，以下代码应用过滤器来查找其邮政编码与地址的邮政编码相匹配且其地址范围包含该地址的街道号码的街道路段。

# Address to be geocoded
address = "2709 Broadway, New York, NY 10025"
zipcode = address.split(" ")[-1]
street_num = address.split(" ")[0] # Find street segments whose left side address ranges contain the street number
potentials = ny_streets.loc[ny_streets["l_low_hn"] street_num]
potentials = potentials.loc[potentials["l_high_hn"] > street_num]
# Find street segments whose zipcode matches the address'
potentials = potentials.loc[potentials["l_zip"] == zipcode]

这会将列表范围缩小到下面所示的一条街段。

XXXXX

最后的任务是确定地址位于该行的位置。这是通过将街道号码放置在段的地址范围内，进行标准化以确定地址沿线的距离，并将该常量应用于线端点的坐标以获得地址的坐标来完成的。以下代码概述了此过程。

import numpy as np
from shapely.geometry import Point # Calculate how far along the street to place the point
denom = ( potentials["l_high_hn"].astype(float) - potentials["l_low_hn"].astype(float)
).values[0]
normalized_street_num = ( float(street_num) - potentials["l_low_hn"].astype(float).values[0]
) / denom # Define a point that far along the street
# Move the line to start at (0,0)
pizza = np.array(potentials["geometry"].values[0].coords[1]) - np.array( potentials["geometry"].values[0].coords[0]
)
# Multiply by normalized street number to get coordinates on line
pizza = pizza * normalized_street_num
# Add starting segment to place line back on the map
pizza = pizza + np.array(potentials["geometry"].values[0].coords[0])
# Convert to geometry array for geopandas
pizza = gpd.GeoDataFrame( {"address": , "geometry": [Point(pizza[0], pizza[1])]}, crs=ny_streets.crs, geometry="geometry",
)

完成地址地理编码后，现在可以在地图上绘制该披萨店的位置以了解其位置。由于上面的代码查看了与街道路段左侧相关的信息，因此实际位置将稍微偏向道路左侧建筑物中绘制点的左侧。你终于知道在哪里可以买到披萨了。

XXXXX

此过程涵盖了最常被称为地理编码的内容，但这并不是该术语的唯一使用方式。您还可能会看到地理编码是指将地标名称转换为坐标、将邮政编码转换为坐标或将坐标转换为 GIS 矢量的过程。您甚至可能听到反向地理编码（稍后将介绍）称为地理编码。包含这些内容的地理编码的更宽松的定义是“位置和地理坐标的近似自然语言描述之间的转换”。因此，每当您需要在这两种数据之间移动时，请考虑将地理编码作为解决方案。

作为每当需要对地址进行地理编码时重复此过程的替代方法，可以使用各种 API 端点，例如美国人口普查局地理编码器和谷歌地理编码 API，免费提供准确的地理编码服务。一些付费选项，例如 Esri 的 ArcGIS, 吉奥科迪奥及 Smarty的甚至为选定的地址提供屋顶精度，这意味着返回的坐标恰好落在建筑物的屋顶上，而不是附近的街道上。以下部分以美国人口普查局地理编码器为例概述了如何使用这些服务将地理编码融入您的数据管道。

为了在地理编码时获得尽可能高的准确性，您应该始终首先确保您的地址格式符合您所选服务的标准。每种服务之间的格式略有不同，但常见的格式是 USPS 格式“PRIMARY# STREET, CITY, STATE, ZIP”，其中 STATE 是缩写代码，PRIMARY# 是街道号码，所有提及的套房号码、建筑物号码和邮政信箱被删除。

地址格式化后，您需要将其提交给 API 进行地理编码。对于美国人口普查局地理编码器，您可以通过单行地址处理选项卡手动提交地址，也可以使用提供REST API 以编程方式提交地址。美国人口普查局地理编码器还允许您使用批量地理编码器对整个文件进行地理编码，并使用基准参数指定数据源。要对之前的披萨店进行地理编码，此链接可用于将地址传递给 REST API，这可以在 Python 中使用以下代码完成。

# Submit the address to the U.S. Census Bureau Geocoder REST API for processing
response = requests.get( "https://geocoding.geo.census.gov/geocoder/locations/onelineaddress?address=2709+Broadway%2C+New+York%2C+NY+10025&benchmark=Public_AR_Current&format=json"
).json()

返回的数据是一个 JSON 文件，可以轻松解码为 Python 字典。它包含一个“tigerLineId”字段，可用于匹配最近街道的 shapefile，一个“side”字段，可用于确定地址位于该街道的哪一侧，以及“fromAddress”和“toAddress”字段其中包含街道段的地址范围。最重要的是，它包含一个“坐标”字段，可用于在地图上定位地址。以下代码从 JSON 文件中提取坐标并将其处理为 GeoDataFrame 以为空间分析做好准备。

# Extract coordinates from the JSON file
coords = response["result"]["addressMatches"][0]["coordinates"]
# Convert coordinates to a Shapely Point
coords = Point(coords["x"], coords["y"])
# Extract matched address
matched_address = response["result"]["addressMatches"][0]["matchedAddress"]
# Create a GeoDataFrame containing the results
pizza_point = gpd.GeoDataFrame( {"address": [matched_address], "geometry": coords}, crs=ny_streets.crs, geometry="geometry",
)

对该点进行可视化显示，它稍微偏离手动地理编码点左侧的道路。

XXXXX

反向地理编码是获取地理坐标并将其与地理区域的自然语言描述进行匹配的过程。如果应用正确，它是在数据科学工具包中附加外部数据的最强大的技术之一。反向地理编码的第一步是确定您的目标地理位置。这是将包含坐标数据的区域。一些常见的例子是人口普查区、邮政编码和城市。第二步是确定该点位于这些区域中的哪一个（如果有）。当使用公共区域时，美国人口普查地理编码器可通过对 REST API 请求进行少量更改来逆向地理编码。链接了确定哪些人口普查地理区域包含之前的披萨店的请求此处。可以使用与之前相同的方法来处理该查询的结果。然而，创造性地定义区域以满足分析需求并手动对其进行反向地理编码开辟了许多可能性。

要手动反转地理编码，您需要确定区域的位置和形状，然后确定该点是否位于该区域的内部。确定一个点是否在多边形内部实际上是一个相当困难的问题，但是光线投射算法，其中从该点开始并沿某个方向无限传播的射线与区域边界相交的次数为奇数次（如果它位于区域内部），否则为偶数次（Shimrat），可用于在大多数情况下求解该问题案例。对于数学爱好者来说，这实际上是若尔当曲线定理（霍施）。请注意，如果您使用来自世界各地的数据，光线投射算法实际上可能会失败，因为光线最终会环绕地球表面并变成一个圆圈。在这种情况下，您必须找到该区域和该点的绕数 (Weisstein)。如果绕数不为零，则该点位于该区域内。幸运的是，Python 的 geopandas 库提供了定义多边形区域内部和测试点是否位于其中所需的功能，而无需所有复杂的数学运算。

虽然手动地理编码对于许多应用程序来说可能过于复杂，但手动反向地理编码可以为您的技能组合提供实用的补充，因为它允许您轻松地将点与高度自定义的区域相匹配。例如，假设您想带一片披萨去公园野餐。您可能想知道披萨店是否位于公园附近。纽约市为其公园提供 shapefile，作为公园属性数据集（纽约公园开放数据团队），也可以使用以下代码通过 SODA API 访问它们。

# Pull NYC park shapefiles
parks = gpd.read_file( BytesIO( requests.get( "https://data.cityofnewyork.us/resource/enfh-gkve.geojson?$limit=5000" ).content )
)
# Limit to parks with green area for a picnic
parks = parks.loc[ parks["typecategory"].isin( [ "Garden", "Nature Area", "Community Park", "Neighborhood Park", "Flagship Park", ] )
]

可以将这些公园添加到可视化中，以查看披萨店附近有哪些公园。

XXXXX

附近显然有一些选项，但使用形状文件和点计算距离可能很困难，而且计算成本很高。相反，可以应用反向地理编码。如上所述，第一步是确定要将点附加到的区域。在本例中，该区域“距离纽约市的一个公园 1/2 英里”。第二步是计算该点是否位于某个区域内，这可以使用前面提到的方法或通过应用 geopandas 中的“contains”函数以数学方式完成。以下代码用于在测试查看哪些公园的缓冲区域现在包含该点之前向公园的边界添加 1/2 英里的缓冲区。

# Project the coordinates from latitude and longitude into meters for distance calculations
buffered_parks = parks.to_crs(epsg=2263)
pizza_point = pizza_point.to_crs(epsg=2263)
# Add a buffer to the regions extending the border by 1/2 mile = 2640 feet
buffered_parks = buffered_parks.buffer(2640)
# Find all parks whose buffered region contains the pizza parlor
pizza_parks = parks.loc[buffered_parks.contains(pizza_point["geometry"].values[0])]

该缓冲区显示附近的公园，在下图中以蓝色突出显示

XXXXX

成功进行反向地理编码后，您了解到披萨店半英里范围内有 8 个公园可供您野餐。享受那片。

XXXXX
j4p4n 制作的披萨片

来源

乔希·劳勒和彼得·希斯。 ESRM 250：森林资源地理信息系统简介。 GIS 的定义，12 年 2009 月 XNUMX 日，西雅图华盛顿大学。课堂讲座。 https://courses.washington.edu/gis250/lessons/introduction_gis/definitions.html
中海集运酒吧。纽约开放数据。 https://data.cityofnewyork.us/City-Government/road/svwp-sbcd
美国人口普查局地理编码器文档。 2022 年 XNUMX 月。 https://geocoding.geo.census.gov/geocoder/Geocoding_Services_API.pdf
Shimrat, M.，“算法 112：点相对于多边形的位置”1962 年， ACM的通讯 第 5 卷第 8 期，1962 年 XNUMX 月。 https://dl.acm.org/doi/10.1145/368637.368653
Hosch，William L.“乔丹曲线定理”。 大英百科全书，13年2018月XNUMX日， https://www.britannica.com/science/Jordan-curve-theorem
Weisstein，Eric W.“轮廓缠绕数”。从 MathWorld–Wolfram 网络资源。 https://mathworld.wolfram.com/ContourWindingNumber.html
纽约市公园开放数据团队。公园物业。 14 年 2023 月 XNUMX 日。 https://nycopendata.socrata.com/Recreation/Parks-Properties/enfh-gkve
j4p4n，“披萨片。” 从 开放美工图库. https://openclipart.org/detail/331718/pizza-slice