数据科学家的地理编码

本文介绍地理编码作为数据科学流程的一部分它涵盖了手动和基于API的地理编码，并提供了一个有趣而引人入胜的示例

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

什么是地理编码？

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

• 与PandasAI交流您的数据集
• 使用Amazon SageMaker在自定义数据集上实现多目标跟踪解决方案
• 使用Amazon SageMaker Python SDK从Amazon SageMaker离线特征存储构建适用于机器学习的数据集

import geopandas as gpd
import requests
from io import BytesIO

# 从SODA API请求数据
req = requests.get(
    "https://data.cityofnewyork.us/resource/gdww-crzy.geojson?stname_lab=BROADWAY"
)
# 将其转换为字节流
reqstrm = BytesIO(req.content)
# 将流读取为GeoDataFrame
ny_streets = gpd.read_file(reqstrm)

此查询的结果超过700个，但这并不意味着您必须检查700个街道才能找到您的比萨。通过可视化数据，您可以看到有3个主要的Broadway街道和一些较小的街道。

这是因为每个街道被分成相应地块的部分，允许更细粒度地查看数据。流程的下一步是确定地址位于这些部分中的哪一个，使用邮政编码和街道号码。数据集中的每个街道段都包含左侧和右侧建筑物地址的地址范围。同样，每个段都包含左侧和右侧街道的邮政编码。为了查找正确的段，以下代码应用过滤器查找其邮政编码与地址的邮政编码相匹配且其地址范围包含地址的街道号码的街道段。

# 要进行地理编码的地址
address = "2709 Broadway, New York, NY 10025"
zipcode = address.split(" ")[-1]
street_num = address.split(" ")[0]

# 找到左侧地址范围包含街道号的街道段
potentials = ny_streets.loc[ny_streets["l_low_hn"] < street_num]
potentials = potentials.loc[potentials["l_high_hn"] > street_num]
# 找到邮政编码与地址匹配的街道段
potentials = potentials.loc[potentials["l_zip"] == zipcode]

这将缩小列表到下面看到的一个街道段。

最后的任务是确定地址在这条线路上的位置。这是通过将街道号码放入线段的地址范围内来完成的，规范化以确定地址应该在线段上的位置，然后将该常数应用于线段端点的坐标以获得地址的坐标。以下代码概述了这个过程。

import numpy as np
from shapely.geometry import 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

# 定义一个点
# 将线段移动到 (0,0) 开始
pizza = np.array(potentials["geometry"].values[0].coords[1]) - np.array(
    potentials["geometry"].values[0].coords[0]
)
# 乘以规范化的街道号码以获取在线上的坐标
pizza = pizza * normalized_street_num
# 添加起始线段以将线段放回地图上
pizza = pizza + np.array(potentials["geometry"].values[0].coords[0])
# 将其转换为 Geopandas 的几何数组
pizza = gpd.GeoDataFrame(
    {"address": [address], "geometry": [Point(pizza[0], pizza[1])]},
    crs=ny_streets.crs,
    geometry="geometry",
)

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

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

作为重复此过程以进行地理编码的替代方法，许多 API 端点，例如美国人口普查局地理编码器和 Google 地理编码 API，提供了免费的准确地理编码服务。一些付费选项，例如 Esri 的 ArcGIS、Geocodio 和 Smarty，甚至为某些地址提供屋顶精度，这意味着返回的坐标恰好落在建筑物的屋顶上，而不是在附近的街道上。以下部分概述了如何使用这些服务将地理编码适配到您的数据管道中，以美国人口普查局地理编码器为例。

如何进行地理编码

为了在进行地理编码时获得最高可能的准确性，您应始终确保地址格式符合所选服务的标准。这将在每个服务之间略有不同，但常见格式是 USPS 格式“PRIMARY# STREET, CITY, STATE, ZIP”，其中 STATE 是缩写代码，PRIMARY# 是街道号码，所有套房号码、建筑物号码和 PO 箱号码都已删除。

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

# 提交地址到美国人口普查局地理编码器 REST API 进行处理
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”字段，可用于匹配最近街道的形状文件，一个“side”字段，可用于确定地址在该街道的哪一侧，以及“fromAddress”和“toAddress”字段，包含街段的地址范围。最重要的是，它包含一个“coordinates”字段，可用于在地图上定位地址。以下代码从JSON文件中提取坐标并将其处理成GeoDataFrame，为空间分析做准备。

# 从JSON文件中提取坐标
coords = response["result"]["addressMatches"][0]["coordinates"]
# 将坐标转换为Shapely点
coords = Point(coords["x"], coords["y"])
# 提取匹配的地址
matched_address = response["result"]["addressMatches"][0]["matchedAddress"]
# 创建包含结果的GeoDataFrame
pizza_point = gpd.GeoDataFrame(
    {"address": [matched_address], "geometry": coords},
    crs=ny_streets.crs,
    geometry="geometry",
)

将这个点可视化后发现它略微偏离手动地理编码的点的左侧。

如何反向地理编码

反向地理编码是将地理坐标与地理区域的自然语言描述相匹配的过程。当正确应用时，它是附加外部数据到数据科学工具箱中最强大的技术之一。反向地理编码的第一步是确定目标地理区域。这是将包含您的坐标数据的区域。一些常见的例子是人口普查区、邮政编码和城市。第二步是确定该点是否在这些区域中的任何一个中。当使用常见区域时，可以通过对REST API请求进行小的更改，使用美国人口普查地理编码器进行反向地理编码。这里链接了一个用于确定哪些人口普查地理区域包含之前的比萨店的请求。此查询的结果可以使用与之前相同的方法处理。但是，创造性地定义地区以适应分析需要并手动反向地理编码可以开启许多可能性。

要手动反向地理编码，您需要确定区域的位置和形状，然后确定点是否在该区域内部。确定一个点是否在一个多边形内部实际上是一个相当困难的问题，但是可以使用射线投射算法，其中从点开始并沿着方向无限地前进的射线与该区域的边界相交奇数次，如果在区域内部，则交叉偶数次（Shimrat），在大多数情况下可以解决。对于数学倾向的人来说，这实际上是Jordan曲线定理（Hosch）的直接应用。请注意，如果您使用来自全球的数据，则射线投射算法实际上可能会失败，因为射线最终将绕过地球表面并变成一个圆。在这种情况下，您将不得不找到该区域和该点的绕数（Weisstein）。如果绕数不为零，则该点在该区域内。幸运的是，Python的geopandas库提供了定义多边形区域内部和测试点是否在其中的功能。

虽然手动地理编码可能对许多应用程序来说过于复杂，但手动反向地理编码可以成为您技能集的实用补充，因为它允许您轻松将您的点与高度定制的区域匹配。例如，假设您想把比萨饼拿到公园里野餐。您可能想知道比萨店是否与公园的距离很近。纽约市提供了其公园的形状文件作为Parks Properties数据集的一部分（NYC Parks Open Data Team），也可以使用以下代码通过其SODA API访问。

# 获取纽约市公园的形状文件
parks = gpd.read_file(
    BytesIO(
        requests.get(
            "https://data.cityofnewyork.us/resource/enfh-gkve.geojson?$limit=5000"
        ).content
    )
)
# 限制公园具有野餐绿地
parks = parks.loc[
    parks["typecategory"].isin(
        [
            "Garden",
            "Nature Area",
            "Community Park",
            "Neighborhood Park",
            "Flagship Park",
        ]
    )
]

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

显然附近有一些选项，但是使用形状文件和点来计算距离可能会很困难且计算量大。相反，可以应用反向地理编码。第一步，如上所述，是确定要将点连接到的区域。在这种情况下，该区域为“纽约市公园半英里范围内”。第二步是计算点是否位于区域内，这可以通过以前提到的方法进行数学计算或通过在geopandas中应用“包含”函数来完成。以下代码用于将1/2英里缓冲区添加到公园的边界上，然后测试哪些公园的缓冲区现在包含该点。

# 将纬度和经度的坐标投影到米上进行距离计算
buffered_parks = parks.to_crs(epsg=2263)
pizza_point = pizza_point.to_crs(epsg=2263)
# 将区域扩展1/2英里= 2640英尺
buffered_parks = buffered_parks.buffer(2640)
# 找出所有缓冲区包含比萨店的公园
pizza_parks = parks.loc[buffered_parks.contains(pizza_point["geometry"].values[0])]

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

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

参考资料

1. Lawler，Josh和Schiess，Peter。 ESRM 250：森林资源中的地理信息系统介绍。 GIS的定义，2009年2月12日，华盛顿大学，西雅图。课堂讲座。 https://courses.washington.edu/gis250/lessons/introduction_gis/definitions.html
2. CSCL PUB。 纽约开放数据。 https://data.cityofnewyork.us/City-Government/road/svwp-sbcd
3. 美国人口普查局地理编码器文档。 2022年8月。 https://geocoding.geo.census.gov/geocoder/Geocoding_Services_API.pdf
4. Shimrat，M.，“算法112：点相对于多边形的位置”1962年，ACM通讯卷5号8月，
   1962. https://dl.acm.org/doi/10.1145/368637.368653
5. Hosch，William L。 “Jordan curve theorem”。 大英百科全书，2018年4月13日，https://www.britannica.com/science/Jordan-curve-theorem
6. Weisstein，Eric W。 “等高线绕数。”来自MathWorld- Wolfram Web资源。 https://mathworld.wolfram.com/ContourWindingNumber.html
7. 纽约市公园开放数据团队。 公园属性。 2023年4月14日。 https://nycopendata.socrata.com/Recreation/Parks-Properties/enfh-gkve
8. j4p4n，“比萨片”。 来自OpenClipArt。 https://openclipart.org/detail/331718/pizza-slice

Evan Miller是Tech Impact的数据科学研究员，他使用数据支持非营利组织和政府机构的社会使命。以前，Evan在中密歇根大学使用机器学习来训练自动驾驶汽车。