データサイエンティストのためのジオコーディング

Geocoding for Data Scientists.

データサイエンティストがデータの「場所」に関するすべてを知る必要がある場合、彼らはしばしば地理情報システム（GIS）に頼ります。GISは、さまざまな目的に役立つ複雑な技術とプログラムのセットですが、ワシントン大学は「地理情報システムは、地球表面の機能についての知識を伝える目的で関連付けられたまたは接続された複雑な配置の集合」というかなり包括的な定義を提供しています（Lawler et al）。GISは、取得から可視化まで、空間データを処理するための広範な技術を包括しており、GIS専門家でなくても有用なツールが多数あります。この記事では、Pythonを使用したジオコーディングの包括的な概要と、いくつかの実用的なアプリケーションのデモを提供します。具体的には、ニューヨーク市、ニューヨーク州にあるピザ店の正確な場所を住所を使用して決定し、近くの公園に関するデータに接続します。デモではPythonコードを使用していますが、コアコンセプトは多くのプログラミング環境に適用でき、ジオコーディングをワークフローに統合するための基盤を提供し、データを空間データに変換し、より複雑な地理的分析のための扉を開きます。

ジオコーディングとは何ですか？

ジオコーディングは、住所データをマッピング座標に変換することを最も一般的に定義されます。通常、これには、住所で通りの名前を検出し、データベース内の現実世界の相当する境界にその通りを一致させ、その通りのどこに住所を配置するかを推定することが含まれます。例として、ニューヨークのブロードウェイにあるピザ店の住所である「2709 Broadway, New York, NY 10025」の単純な手動ジオコードのプロセスを実行してみましょう。最初のタスクは、住所の場所の道路システムに適したシェープファイルを見つけることです。この場合、住所の市と州は「ニューヨーク、ニューヨーク州」となります。幸いなことに、ニューヨーク市は詳細な道路情報をNYC Open Dataページ（CSCL PUB）に公開しています。次に、「ブロードウェイ」という通り名を調べます。これで、ニューヨーク市内の「ブロードウェイ」という名前の通りに住所がある可能性があることがわかります。以下のPythonコードを実行して、NYC Open Data SODA APIから「Broadway」という名前のすべての通りをクエリできます。

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つの主要なブロードウェイ通りといくつかの小さな通りがあることがわかります。

これは、各通りがおおよそ1ブロックに対応するセクションに分割されているためで、より細かいデータを見ることができます。プロセスの次のステップは、ZIPコードとストリート番号を使用して、住所がどのセクションにあるかを正確に特定することです。データセットの各通りセグメントには、通りの左側と右側の建物の住所範囲が含まれています。同様に、各セグメントには、通りの左側と右側のZIPコードが含まれています。正しいセグメントを特定するには、以下のコードを適用して、住所のZIPコードに一致する通りセグメントを見つけ、住所のストリート番号を含むアドレス範囲を持つ通りセグメントを見つけます。

# ジオコードする住所
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]
# 住所のZIPコードに一致する通りセグメントを見つける
potentials = potentials.loc[potentials["l_zip"] == zipcode]

これにより、下の1つの通りセグメントに絞り込まれます。

最終的なタスクは、住所がこのライン上のどこにあるかを決定することです。これは、セグメントの住所範囲にストリート番号を配置し、正規化して、住所がライン上のどの程度の位置にあるかを決定し、その定数をラインの端点の座標に適用して住所の座標を取得することで行われます。次のコードはこのプロセスを概説しています。

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": [address], "geometry": [Point(pizza[0], pizza[1])]},
    crs=ny_streets.crs,
    geometry="geometry",
)

住所のジオコーディングが完了したので、このピザ店の場所を地図上にプロットして、その場所を把握することができます。上記のコードは、通りの左側に関する情報を見たため、実際の場所は道路の左側の建物内に、プロットされた点よりもやや左になります。これで、ピザを手に入れることができる場所が分かりました。

このプロセスは、一般的にジオコーディングと呼ばれるものですが、これが用語で使われる唯一の方法ではありません。ランドマーク名を座標に変換したり、郵便番号を座標に変換したり、座標をGISベクトルに変換したりするプロセスも、ジオコーディングとして言及されることがあります。さらに、後に説明するリバースジオコーディング（逆ジオコーディング）もジオコーディングとして言及されることがあります。場所のおおよその自然言語表現と地理座標の間を移動する必要がある場合は、ジオコーディングを解決策として検討してください。

住所のジオコーディングを繰り返す代わりに、米国国勢調査局ジオコーダーやGoogle Geocoding APIなどのAPIエンドポイントなど、正確なジオコーディングサービスを無償で提供するものがあります。EsriのArcGIS、Geocodio、Smartyなどの有料オプションでは、選択した住所に対して屋上の正確性を提供する場合があります。これは、返される座標が近くの通りではなく、建物の屋根に正確に着地することを意味します。次のセクションでは、米国国勢調査局ジオコーダーを例にして、これらのサービスを使用してジオコーディングをデータパイプラインに適合させる方法について説明します。

ジオコーディングの方法

ジオコーディング時に最高の精度を得るためには、常に選択したサービスの標準に合わせて住所のフォーマットを確認する必要があります。これは、各サービスごとにわずかに異なりますが、一般的なフォーマットは「PRIMARY＃STREET、CITY、STATE、ZIP」というUSPSフォーマットです。ここで、STATEは略語コード、PRIMARY＃は通りの番号で、スイート番号、建物番号、POボックスのすべての言及が削除されます。

住所をフォーマットしたら、APIに住所をジオコードするために送信する必要があります。米国国勢調査局ジオコーダーの場合、One Line Address Processingタブを介して住所を手動で送信するか、提供された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」というフィールドが含まれており、最も近いストリートのシェープファイルに一致させるために使用できます。また、「side」というフィールドが含まれており、そのアドレスがそのストリートのどちら側にあるかを決定するために使用できます。そして、ストリートセグメントの住所範囲を含む「fromAddress」と「toAddress」というフィールドが含まれています。最も重要なのは、「coordinates」というフィールドであり、これを使用して地図上のアドレスを検出できます。次のコードは、JSONファイルから座標を抽出し、ジオデータフレームに処理して、空間分析の準備をします。

# JSONファイルから座標を抽出する
coords = response["result"]["addressMatches"][0]["coordinates"]
# 座標をShapely Pointに変換する
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",
)

このポイントを可視化すると、手動でジオコーディングされたポイントよりも道路から左に少し外れていることがわかります。

逆ジオコーディングの方法

逆ジオコーディングとは、地理座標を自然言語の地理的範囲の説明に一致させるプロセスです。正しく適用されると、データサイエンスのツールキットで最も強力なテクニックの1つです。逆ジオコーディングの最初のステップは、ターゲット地理情報を決定することです。これは、座標データを含む領域です。一般的な例には、国勢調査トラクト、郵便番号、都市などがあります。2番目のステップは、その地点がどの地域に含まれるかを決定することです。一般的な地域を使用する場合、米国国勢調査ジオコーダーを使用してREST APIリクエストを微調整することで、逆ジオコーディングを実行できます。以前のピザ店が含まれる国勢調査地理情報を決定するためのリクエストは、こちらにリンクされています。このクエリの結果は、前と同じ方法で処理できます。ただし、解析ニーズに合わせて領域を創造的に定義し、手動で逆ジオコーディングを行うことで、多くの可能性が開けます。

手動で逆ジオコーディングを行うには、領域の場所と形状を決定し、その点がその領域の内部にあるかどうかを判断する必要があります。点が多角形の内部にあるかどうかを判断することは、実際にはかなり難しい問題ですが、レイキャスティングアルゴリズム（Shimrat）を使用することで、ほとんどの場合解決できます。レイキャスティングアルゴリズムとは、点から始まり、方向を無限に伸ばし、その方向の直線が領域の境界線と奇数回交差する場合は領域の内部にあり、偶数回交差する場合は領域の外部にあるというアルゴリズムです（Hosch）。数学的に傾斜しやすい人にとっては、これは実際にはJordanの曲線定理（Hosch）の直接的な応用です。なお、世界中のデータを使用する場合、レイキャスティングアルゴリズムは実際には失敗することがあります。これは、レイが最終的に地球の表面を包み込み、円になるためです。この場合、代わりに領域と点の巻き数（Weisstein）を見つける必要があります。巻き数がゼロでない場合、点は領域の内部にあります。幸いなことに、Pythonのgeopandasライブラリは、複雑な数学を必要とせずに多角形領域の内部を定義し、点がそれに含まれているかどうかをテストするために必要な機能を提供しています。

手動ジオコーディングは多くのアプリケーションにとって複雑すぎる場合がありますが、手動逆ジオコーディングは、高度にカスタマイズされた領域に簡単にポイントを一致させることができるため、スキルセットに実用的な追加となります。たとえば、ピザを持って公園に行ってピクニックをしたいと思った場合、ピザ店が公園から近いかどうかを知りたい場合があります。ニューヨーク市は、Parks Propertiesデータセットの一部として、公園のシェープファイルを提供しており、次のコードを使用してSODA APIを介してアクセスすることもできます。

# NYC公園のシェープファイルを取得する
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",
        ]
    )
]

これらの公園は、視覚化に追加され、ピザ店の近くにどのような公園があるかを確認することができます。

明らかにいくつかのオプションが近くにありますが、シェイプファイルとポイントを使用して距離を算出することは困難で、計算コストも高くなります。代わりに、逆ジオコーディングを適用することができます。まず上記で言及したように、ポイントをアタッチするリージョンを決定することが最初のステップです。この場合、リージョンは「ニューヨーク市の公園から半マイルの距離」となります。2番目のステップは、ポイントが領域内にあるかどうかを計算することであり、これは以前に言及した方法を数学的に使用するか、geopandasの「contains」関数を適用することで行うことができます。以下のコードは、テストする前に公園の境界に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

参考文献

Lawler, Josh and Schiess, Peter. ESRM 250: Introduction to Geographic Information Systems in Forest Resources. Definitions of GIS, 12 Feb. 2009, University of Washington, Seattle. Class Lecture. https://courses.washington.edu/gis250/lessons/introduction_gis/definitions.html
CSCL PUB. New York OpenData. https://data.cityofnewyork.us/City-Government/road/svwp-sbcd
U.S. Census Bureau Geocoder Documentation. August 2022. https://geocoding.geo.census.gov/geocoder/Geocoding_Services_API.pdf
Shimrat, M., “Algorithm 112: Position of point relative to polygon” 1962, Communications of the ACM Volume 5 Issue 8, Aug.
1. https://dl.acm.org/doi/10.1145/368637.368653
Hosch, William L.. “Jordan curve theorem”. Encyclopedia Britannica , 13 Apr. 2018, https://www.britannica.com/science/Jordan-curve-theorem
Weisstein, Eric W. “Contour Winding Number.” From MathWorld –A Wolfram Web Resource. https://mathworld.wolfram.com/ContourWindingNumber.html
NYC Parks Open Data Team. Parks Properties. April 14, 2023. https://nycopendata.socrata.com/Recreation/Parks-Properties/enfh-gkve
j4p4n, “Pizza Slice.” From OpenClipArt . https://openclipart.org/detail/331718/pizza-slice

Evan Millerは、社会貢献を目的とした非営利および政府機関を支援するためにデータを使用するTech Impactのデータサイエンスフェローです。以前、EvanはCentral Michigan Universityで自律型車両のトレーニングに機械学習を使用しました。

We will continue to update VoAGI; if you have any questions or suggestions, please contact us!

Data science

Was this article helpful?

93 out of 132 found this helpful