library(tidyverse)
library(sf)
library(tmap)
library(tmaptools)
library(tidygeocoder)ジオコーディング
ジオコーディング
ジオコーディングとは
住所を緯度経度に変換することは可能でしょうか? 実際に地域分析を行う場面においては、位置データとして住所情報のみが得られる場合が多く、 これを分析する際に、住所を緯度経度に変換する必要に迫られる事態に頻繁に遭遇します。
この、住所を緯度経度に変換することを、ジオコーディング(Geocoding)と呼びます。 そして緯度経度を住所に変換することは、逆ジオコーディング(Reverse Geocoding)と呼びます。
Rのみでのジオコーディングはかなり煩雑なので, ここでは,外部のサーピスやアプリを使ったいくつかのジオコーディングの方法を紹介します。
2023年6月頃、一部界隈で日本の住所がヤバいという話題が盛り上がりました。 詳しくはリンク先の記事を読んでいただきたいのですが、ざっくりまとめると、日本の住所には表記の揺れやローカルルールなどの問題があり、全国の住所を正規化することが困難だという内容です。 つまり収集した住所データをコンピュータの自動処理にかける前の手間がかかることを意味するのですが、このことは最近になって判明したことではなく以前からよく知られたことでした。 しかし2023年6月のTwitterでの発言をきっかけに盛り上がりました。 というわけで、ジオコーディングをサービスとして行う民間企業も以前からあったのですが、住所正規化に商機を見出したのか、この頃から各社とも積極的にサービスをPRしているような気がします。
ただしこれらのサービスの多くは、主に法人を対象としたサービスとなっており、個人として利用することのハードルは、費用も含めて、かなり高いです。
外部のサービスとのデータのやり取りに、API1を使う方法と、ファイルを使う方法があります。
Rから外部のAPIを利用するパッケージを利用することで、簡単にジオコーディングを行うことができます。
また、外部サービスによるジオコードの結果をRに読み込み, 空間データとして分析に使用する方法も併せて紹介します。
APIを利用する方法では、tidygeocoderパッケージを利用する方法を解説します。 ファイルを利用する方法では、国土地理院のWeb地図「地理院地図」を使ったジオコーディングの方法を紹介します。
使用するデータ
今回の演習では、佐賀県内にある「スターバックスコーヒー」の店舗の住所を使います。 はじめに、以下のCSVファイルをダウンロード(右クリックから[名前をつけてリンクを保存…]を選択)して、ダウンロードしたファイルをプロジェクトのdataフォルダに移動してください。
これは、佐賀県内のスターバックス店舗の住所一覧を、公式ウェブサイトの情報をもとに作成したものです。 ファイルの中身は、14店舗の店名と住所だけです。

tidygeocoderによるジオコーディング
tidygeocoderとは
公式サイトには、以下のように説明されています。
Tidygeocoder makes getting data from geocoding services easy. A unified high-level interface is provided for a selection of supported geocoding services and results are returned in tibble (dataframe) format.
tidygeocoderよるジオコーディング
なにはともあれ、tidygeocoderパッケージをロードしましょう。
使い方は極めて簡単です。geo関数に住所の文字列を渡すだけです。 method引数で使用するAPIを指定するのですが、ここではarcgisを使っています。
geo('佐賀県佐賀市本庄町1', method = 'arcgis')Passing 1 address to the ArcGIS single address geocoder
Query completed in: 0.6 seconds
# A tibble: 1 × 3
address lat long
<chr> <dbl> <dbl>
1 佐賀県佐賀市本庄町1 33.2 130.
tidygeocoderが対応しているジオコーディングAPIの一覧は、公式サイトをご覧ください。 多くのサービスではAPI Keyが必要となるのですが、API Keyの取得にはユーザー登録が必要となり手続きが煩雑なため、ここではAPI Keyの必要ないサービスであるArcGISを利用しています。
geo関数は、データフレームを返します。 address列に加えて、lat列(latitude:緯度)、long列(longitude:経度)のデータからなるデータフレームになっています。
住所データとしては、1つの文字列だけでなく、2つ以上の文字列からなるベクトルを渡すこともできます。
address <-
c("佐賀県佐賀市本庄町1",
"佐賀県佐賀市鍋島5丁目1-1",
"佐賀県西松浦郡有田町大野乙2441-1")
geo_addr <- geo(address, method = 'arcgis')Passing 3 addresses to the ArcGIS single address geocoder
Query completed in: 1.9 seconds
geo_addr# A tibble: 3 × 3
address lat long
<chr> <dbl> <dbl>
1 佐賀県佐賀市本庄町1 33.2 130.
2 佐賀県佐賀市鍋島5丁目1-1 33.3 130.
3 佐賀県西松浦郡有田町大野乙2441-1 33.2 130.
このように、geo関数は、入力した住所のベクトルの長さと同じ行数を持つデータフレームを返します。
ただしこのままでは、ただの経度と緯度の数値データからなるデータフレームなので、これを空間データに変換する必要があります。
geo_addr <-
geo_addr |>
st_as_sf(coords = c("long", "lat"), crs = 4326)
geo_addrSimple feature collection with 3 features and 1 field
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 129.8872 ymin: 33.17706 xmax: 130.2913 ymax: 33.28433
Geodetic CRS: WGS 84
# A tibble: 3 × 2
address geometry
* <chr> <POINT [°]>
1 佐賀県佐賀市本庄町1 (130.2913 33.24219)
2 佐賀県佐賀市鍋島5丁目1-1 (130.2671 33.28433)
3 佐賀県西松浦郡有田町大野乙2441-1 (129.8872 33.17706)
sf::st_as_sf関数でデータフレームをsfのポイントデータに変換します。 coords引数で、データフレームのどの列を座標データに利用するかを指定し、 CRSはWGS84(EPSGコード:4326)を指定しています。
これを地図に表示してみましょう。
gsi_tile <- "https://cyberjapandata.gsi.go.jp/xyz/pale/{z}/{x}/{y}.png"
tm_shape(geo_addr, bbox = bb(geo_addr, ext = 1.5)) +
tm_symbols() + tm_basemap(gsi_tile) +
tm_crs(6670) + tm_compass() + tm_scalebar()[basemaps] Tiles from "id_1",
"https://cyberjapandata.gsi.go.jp/xyz/pale/{z}/{x}/{y}.png", "c(\"a\", \"b\",
\"c\")", and "" will be projected so details (e.g. text) could appear blurry
This message is displayed once per session.

とても簡単ですね。
データフレームのバッチ処理によるジオコーディング
住所の文字列や、そのベクトルを用意すれば、ジオコードが簡単にできることがわかりました。 tidygeocoderパッケージを使えば、データフレームに対して一括してジオコードを行うこともできます。 先ほどダウンロードしたStarbucks_saga.csvをread_csv関数を使って読み込みましょう。
starbucks <- read_csv('data/StarbucksSaga.csv')Rows: 14 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): 名前, 住所
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
starbucks# A tibble: 14 × 2
名前 住所
<chr> <chr>
1 佐賀南バイパス店 佐賀県佐賀市本庄町袋306-6
2 佐賀大学通り店 佐賀県佐賀市与賀町70-1
3 鳥栖プレミアム・アウトレット店 佐賀県鳥栖市弥生が丘8-1
4 TSUTAYA鳥栖店 佐賀県鳥栖市本鳥栖町537-1
5 鳥栖蔵上町店 佐賀県鳥栖市蔵上町662-5
6 基山パーキングエリア(上り線)店 佐賀県三養基郡基山町小倉2097-1
7 吉野ヶ里歴史公園店 佐賀県神埼市神埼町鶴612-1
8 蔦屋書店武雄市図書館店 佐賀県武雄市武雄町武雄5304-1
9 佐賀武雄店 佐賀県武雄市武雄町昭和277
10 唐津店 佐賀県唐津市和多田西山1-46-1
11 TSUTAYABOOKSTORE伊万里店 佐賀県伊万里市二里町八谷搦1274
12 金立サービスエリア(下り線)店 佐賀県佐賀市金立町金立1197-57
13 佐賀兵庫店 佐賀県佐賀市兵庫町藤木1192-1
14 ゆめタウン佐賀店 佐賀県佐賀市兵庫北5-14-1
店舗名(列名:名前)と住所の2列からなるデータフレームで、佐賀県内のR nrow(starbucks)店舗のデータが入っていることがわかります。 データフレームに対してジオコーディングを行うにはgeocode関数を使います。
starbucks_tg <- starbucks |> geocode(住所, method = 'arcgis')Passing 14 addresses to the ArcGIS single address geocoder
Query completed in: 10.2 seconds
starbucks_tg# A tibble: 14 × 4
名前 住所 lat long
<chr> <chr> <dbl> <dbl>
1 佐賀南バイパス店 佐賀県佐賀市本庄町袋306-6 33.2 130.
2 佐賀大学通り店 佐賀県佐賀市与賀町70-1 33.2 130.
3 鳥栖プレミアム・アウトレット店 佐賀県鳥栖市弥生が丘8-1 33.4 131.
4 TSUTAYA鳥栖店 佐賀県鳥栖市本鳥栖町537-1 33.4 131.
5 鳥栖蔵上町店 佐賀県鳥栖市蔵上町662-5 33.4 130.
6 基山パーキングエリア(上り線)店 佐賀県三養基郡基山町小倉2097-1 33.4 131.
7 吉野ヶ里歴史公園店 佐賀県神埼市神埼町鶴612-1 33.3 130.
8 蔦屋書店武雄市図書館店 佐賀県武雄市武雄町武雄5304-1 33.2 130.
9 佐賀武雄店 佐賀県武雄市武雄町昭和277 33.2 130.
10 唐津店 佐賀県唐津市和多田西山1-46-1 33.4 130.
11 TSUTAYABOOKSTORE伊万里店 佐賀県伊万里市二里町八谷搦1274 33.3 130.
12 金立サービスエリア(下り線)店 佐賀県佐賀市金立町金立1197-57 33.3 130.
13 佐賀兵庫店 佐賀県佐賀市兵庫町藤木1192-1 33.3 130.
14 ゆめタウン佐賀店 佐賀県佐賀市兵庫北5-14-1 33.3 130.
geocode関数の引数には、入力のデータフレーム(この例ではstarbucks)の住所データが入っている列名(この場合は住所)を指定します。 出力されたデータフレームであるstarbuks_tgには、経度(lat)と緯度(long)のデータ列が追加されていることが確認できると思います。
これまでと同じように、これを空間データ(sf)に変換して地図に描画してみましょう。
starbucks_tg <- starbucks_tg |>
st_as_sf(coords = c("long", "lat"), crs = 4326)
tm_shape(starbucks_tg, bbox = bb(starbucks_tg, ext = 1.3)) +
tm_symbols() + tm_basemap(gsi_tile) +
tm_crs(6670) + tm_compass() + tm_scalebar()
地理院地図によるジオコーディング
地理院地図とは
地理院地図は、国土地理院が提供する、地図や空中写真などが閲覧できるWebサービスです。 公式サイトには、以下のように説明されています。
地理院地図とは、地形図、写真、標高、地形分類、災害情報など、国土地理院が捉えた日本の国土の様子を発信するウェブ地図です。3Dで見ることもできます。また、地形断面図の作成や新旧の写真を比較する機能なども備えています。
地理院地図によるジオコーディング
地理院地図をブラウザで開いてください。 次のような画面が表示されると思います。

このブラウザのウィンドウに、CSVファイル(Starbucks_saga.csv)をドラッグ&ドロップします(ブラウザにファイルをドラッグ&ドロップするという経験はあまりないかもしれません。最初は混乱するかもしれませんが、大丈夫です)。

すると、「CSVファイル読込」というウィンドウが開きますので、プルダウンメニューから[住所]を選択し、 [上記の内容で読込開始]ボタンをクリックします。

すると、佐賀県内のスターバックスの位置が、赤い丸のアイコン🔴で表示されたのがわかると思います。

住所データが、地図上に表示されたということは、その場所の緯度経度のデータがわかった、つまりジオコードされたことを意味します。
それでは、ジオコードされた緯度経度データをRに取り込む方法について説明します。 地理院地図では、ジオコードされたデータを、KMLファイルまたはGeoJSONファイルとして保存することが可能です。 ここでは、GeoJSONファイルとして保存する方法を説明します。
画面右上の[ツール]をクリックすると、画面右側にメニューが表示されますので、[作図・ファイル]をクリックします。

「作図・ファイル」というウィンドウが開きますので、[Starbucks_saga.csv]にチェックがついていることを確認したら、保存ボタン💾をクリックします。

ファイル形式を選択するウィンドウが開きますので、[GeoJSON形式]の左のラジオボタンを選択し、 下の[上記の内容で保存]をクリックします。

ファイル名を「StarbucksSaga.geojson」にして、[保存]ボタンをクリックします。

保存したGeoJSONファイルをプロジェクトのdataフォルダに移動しましょう。
結果の表示
地理院地図によるジオコーディングを行わなかった場合は、以下のリンクからGeoJSONファイルをダウンロードし、dataフォルダに移動してください。
GeoJSONファイルを、st_read関数で読み込みます。
starbucks_gsi <-
st_read("data/StarbucksSaga.geojson")Reading layer `StarbucksSaga' from data source
`/Users/kzktmr/Documents/SagaU/website/user/geospatial/data/StarbucksSaga.geojson'
using driver `GeoJSON'
Simple feature collection with 14 features and 6 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 129.8677 ymin: 33.1889 xmax: 130.5298 ymax: 33.44066
Geodetic CRS: WGS 84
地理院地図で保存したGeoJSONには、地理院地図で使うためのデータも入っています。
glimpse(starbucks_gsi)Rows: 14
Columns: 7
$ name <chr> "佐賀南バイパス店", "佐賀大学通り店", "鳥栖プレミアム・アウトレット店", "TSUTAYA鳥栖店", …
$ 住所 <chr> "佐賀県佐賀市本庄町袋306-6", "佐賀県佐賀市与賀町70-1", "佐賀県鳥栖市弥生が丘8-1", "佐賀…
$ X_markerType <chr> "Icon", "Icon", "Icon", "Icon", "Icon", "Icon", "Icon", …
$ X_iconUrl <chr> "https://maps.gsi.go.jp/portal/sys/v4/symbols/080.png", …
$ `_iconSize` <list> <20, 20>, <20, 20>, <20, 20>, <20, 20>, <20, 20>, <20, …
$ `_iconAnchor` <list> <10, 10>, <10, 10>, <10, 10>, <10, 10>, <10, 10>, <10, …
$ geometry <POINT [°]> POINT (130.299 33.23963), POINT (130.2923 33.24775…
これはRでは使用しませんので、不要なデータを取り除きましょう。 ここでは、select関数を使ってみます。
starbucks_gsi <- starbucks_gsi |>
select(name, address = 住所, geometry)
glimpse(starbucks_gsi)Rows: 14
Columns: 3
$ name <chr> "佐賀南バイパス店", "佐賀大学通り店", "鳥栖プレミアム・アウトレット店", "TSUTAYA鳥栖店", "鳥栖蔵上…
$ address <chr> "佐賀県佐賀市本庄町袋306-6", "佐賀県佐賀市与賀町70-1", "佐賀県鳥栖市弥生が丘8-1", "佐賀県鳥栖市本…
$ geometry <POINT [°]> POINT (130.299 33.23963), POINT (130.2923 33.24775), PO…
それでは、地図に表示してみましょう。
tm_shape(starbucks_gsi, bbox = bb(starbucks_gsi, ext = 1.3)) +
tm_symbols() + tm_basemap(gsi_tile) +
tm_crs(6670) + tm_compass() + tm_scalebar()
ジオコーディング結果の比較
今回の演習で使用した2つのジオコーディング方法の結果を比較してみましょう。 2つの方法で作ったデータを地図に重ねて表示してみます。 ここではmapviewを使いましょう。
tm_shape(starbucks_tg, bbox = bb(starbucks_tg, ext = 1.3)) + tm_symbols(fill = "red") +
tm_shape(starbucks_gsi) + tm_symbols(fill = "green") +
tm_basemap(gsi_tile) +
tm_crs(6670) + tm_compass() + tm_scalebar()
これをみると、「金立サービスエリア(下り線)店」と「吉野ヶ里歴史公園店」の場所を除けば、2つのジオコーディング結果に大きな違いはないことがわかると思います。金立サービスエリア(下り線)店の位置については、tinygeocoderを使った結果が正確のようです。 また吉野ヶ里歴史公園店の位置については、広大な敷地の吉野ヶ里歴史公園内の店舗ということもあってか、実際の位置から数百メートルずれてしまっています。
郵便番号のジオコーディング
アンケート調査などで、居住地や勤務先の住所を尋ねるのではなく、郵便番号を尋ねることがあります。これは、有効回答率を高めるために、住所よりも回答しやすいと思われる郵便番号を尋ねることで、回答の心理的な障壁を下げるためです。
さて、郵便番号のデータがあるときに、そのジオコーディングはどのようにすればよいでしょうか。
これもtidygeocoderを使うことができます。 method引数にosm(Nominatim:OpenStreetMapのデータを用いたオープンソースのジオコーディングツール)を設定すれば、郵便番号によるジオコーディングが可能です。
postalcode <- c('8498501', '8408502')
postalcode_geo <- geo(postal = postalcode, method = 'osm')Passing 2 addresses to the Nominatim single address geocoder
Query completed in: 2 seconds
postalcode_geo# A tibble: 2 × 3
postalcode lat long
<chr> <dbl> <dbl>
1 8498501 33.3 130.
2 8408502 33.2 130.
postalcode_geo <-
postalcode_geo |>
st_as_sf(coords = c('long', 'lat'), crs = 4326)
tm_shape(postalcode_geo, bbox = bb(postalcode_geo, ext = 1.5)) +
tm_symbols() +
tm_basemap(gsi_tile) +
tm_crs(6670) + tm_compass() + tm_scalebar()
脚注
Application Programming Interface。ソフトウェアとWebサービスとを繋ぐインターフェイスのこと。↩︎