library(tidyverse)
library(sf)
library(mapview)6: ジオコーディング
空間結合
GISの重要な特徴は、複数の空間データを重ねて（オーバーレイと言ったりします）分析できることです。
ここでは、ポイントデータとポリゴンデータの空間結合の例を示します。
ポイントデータにポリゴンデータを空間結合することで、それぞれのポイントデータに対して、そのポイントデータがどのポリゴン領域に含まれるかという情報を付与することができます。 例えば、ある緯度経度が与えられたときに、その緯度経度がどの行政区域に含まれているか調べることができます。
データの準備
はじめに、今日の演習で使用するパッケージをロードします。
ここでは、佐賀県の市町ポリゴンデータと、乱数によって発生させた100個の地点との包含関係を分析します。 まずはじめに、およそ佐賀県に相当する区域内にランダムな100個の点があるような、ポイントデータを作成します。
乱数を発生させる緯度と経度の範囲を知るために、佐賀県の市町ポリゴンデータのGeoJSONファイルを変数sagaに読み込み、 st_bbox関数で、sagaのバウンディングボックス（オブジェクトを包含する最小の矩形領域）を調べます。
佐賀県のGeoJSONファイルは、佐賀県オープンデータカタログサイトからダウンロードし、dataフォルダに入れてください（前回の講義で利用したものと同じなので、前回利用したファイルをコピーしても構いません）。
sf::st_read関数を使って読み込みます。
saga <- st_read("data/410004saga.geojson")Reading layer `410004saga' from data source
`/Users/kzktmr/Documents/LectureNotes/data/410004saga.geojson'
using driver `GeoJSON'
Simple feature collection with 20 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 129.7368 ymin: 32.95054 xmax: 130.5424 ymax: 33.6189
Geodetic CRS: WGS 84バウンディングボックスはsf::st_bbox関数で調べることができます。
saga_bbox <- st_bbox(saga)
saga_bbox xmin ymin xmax ymax
129.73681 32.95054 130.54235 33.61890 佐賀県の経度は129.7368〜130.5424、緯度は32.9505〜33.6189の範囲であることがわかりましたので、この範囲で、 乱数（runif関数）を使って、100個の点からなるポイントデータを作成します。
n <- 100
random_points <-
tibble(x = runif(n, min = saga_bbox$xmin, max = saga_bbox$xmax),
y = runif(n, min = saga_bbox$ymin, max = saga_bbox$ymax))xとyそれぞれに100個の乱数を発生させ、そのxとyを列にもつデータフレームを作成します。 xとyの値の上下限に、さきほど調べたバウンディングボックスを使っています。 どのようなデータなのか確認するために、ggplotを使って、xとyの散布図を描いてみましょう。
random_points |> ggplot(aes(x = x, y = y)) +
geom_point()
このままでは、ただの数値データからなるデータフレームなので、空間データに変換する必要があります。
random_points <- random_points |>
st_as_sf(coords = c("x", "y"), crs = 4326)sf::st_as_sf関数でデータフレームをsfのポイントデータに変換します。 coords引数で、データフレームのどの列を座標データに利用するかを指定し、 CRSはWGS84（EPSGコード：4326）を指定しています。
random_pointsSimple feature collection with 100 features and 0 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 129.7432 ymin: 32.95105 xmax: 130.5284 ymax: 33.61656
Geodetic CRS: WGS 84
# A tibble: 100 × 1
geometry
* <POINT [°]>
1 (129.9084 33.44301)
2 (129.7776 33.45922)
3 (129.9373 33.6145)
4 (130.4801 33.08488)
5 (130.0143 33.16909)
6 (130.1443 32.96391)
7 (130.1597 33.03113)
8 (130.2464 33.24463)
9 (130.2018 33.40191)
10 (130.1549 33.11526)
# ℹ 90 more rowsこれで、100個の点からなるポイントデータrandom_pointsが作成できました。 佐賀県の地図と合わせて、mapview関数で表示してみましょう。
mapview(random_points) + mapview(saga)Warning in cbind(`Feature ID` = fid, mat): number of rows of result is not a
multiple of vector length (arg 1)st_join関数による空間結合
さて、これらポイントデータがどのポリゴンに含まれるかを知りたいわけですが、そのために、ポイントデータにポリゴンデータを結合します。 空間データどうしの結合には、left_join関数ではなく、sfパッケージのst_join関数を使います。
random_points_joined <- st_join(random_points, saga)
glimpse(random_points_joined)Rows: 100
Columns: 5
$ geometry <POINT [°]> POINT (129.9084 33.44301), POINT (129.7776 33.45922), P…
$ KEN_NAME <chr> "佐賀県", NA, NA, NA, "佐賀県", NA, "佐賀県", "佐賀県", "佐賀…
$ GST_NAME <chr> "唐津市", NA, NA, NA, "武雄市", NA, "藤津郡", "佐賀市", "佐賀…
$ CSS_NAME <chr> NA, NA, NA, NA, NA, NA, "太良町", NA, NA, NA, NA, NA, NA, NA,…
$ KEY_CODE <chr> "412023", NA, NA, NA, "412066", NA, "414417", "412015", "4120…これをみると、100個の点それぞれに、市町名や自治体コードが結合されていることがわかります。 ただし、佐賀県の領域の外にある点にはポリゴンの情報が結合されておらず、自治体コードがNAになっています。
そこで、佐賀県の外にある点を、filter関数を使って取り除きましょう。
ヒント
NAは、Not Availableの略で、欠損値を意味します。 値が欠損値であるかどうかを判定するには、is.na関数を使います。 その前に否定の論理演算子である!をつけると、!is.na(x)で「xが欠損値ではない」という意味になります。
random_points_joined <-
st_join(random_points, saga) |>
filter(!is.na(KEY_CODE))
glimpse(random_points_joined)Rows: 36
Columns: 5
$ geometry <POINT [°]> POINT (129.9084 33.44301), POINT (130.0143 33.16909), P…
$ KEN_NAME <chr> "佐賀県", "佐賀県", "佐賀県", "佐賀県", "佐賀県", "佐賀県", "…
$ GST_NAME <chr> "唐津市", "武雄市", "藤津郡", "佐賀市", "佐賀市", "唐津市", "…
$ CSS_NAME <chr> NA, NA, "太良町", NA, NA, NA, NA, NA, NA, NA, NA, NA, "吉野ヶ…
$ KEY_CODE <chr> "412023", "412066", "414417", "412015", "412015", "412023", "…地図に表示して、佐賀県に含まれない点が取り除かれているかどうか確認しましょう。
mapview(random_points_joined) + mapview(saga)ポイントデータとポリゴンデータの両方をそれぞれ同じ配色で描画すると、ポイントデータにポリゴンデータのKEY_CODEの情報が結合されていることが確認できます。
keycode <- sort(saga$KEY_CODE)
saga <- saga |> arrange(KEY_CODE) |>
mutate(KEY_CODE = factor(KEY_CODE, levels = keycode))
random_points_joined <- random_points_joined |>
mutate(KEY_CODE = factor(KEY_CODE, levels = keycode))
mapview(random_points_joined, zcol = "KEY_CODE",
col.regions = rainbow(20)) +
mapview(saga, zcol = "KEY_CODE",
col.regions = rainbow(20))今ここで行った操作は、ポイントデータに市町のデータを結合するものだったのですが、これは言い換えると、ポイントデータに「佐賀県佐賀市」といった住所を紐づけたことに他なりません。 つまり、緯度経度のデータに住所という属性データを結合したというわけです。
今回は市町のポリゴンデータを使いましたが、これをもっと細かいポリゴンデータ、例えば小地域データ（何丁目あるいは大字）や街区レベル（何番何号）のポリゴンデータを使えば、緯度経度を住所に変換できることがわかると思います。
ジオコーディング
ジオコーディングとは
繰り返しになりますが、先ほどの空間結合の例は、緯度経度の情報を住所に変換する基本的な考え方の説明になっています。
それでは、その逆に、住所を緯度経度に変換することは可能でしょうか？ 答えは、住所と緯度経度を含む空間データがあれば、可能だといえます。 実際に地域分析を行う場面においては、位置データとして住所情報のみが得られる場合が多く、 これを分析する際に、住所を緯度経度に変換する必要に迫られる事態に頻繁に遭遇します。
この、住所を緯度経度に変換することを、ジオコーディング（Geocoding）と呼びます。 そして緯度経度を住所に変換することは、逆ジオコーディング（Reverse Geocoding）と呼びます。
Rのみでのジオコーディングはかなり煩雑なので， ここでは，外部のサーピスやアプリを使ったいくつかのジオコーディングの方法を紹介します。
ノート
2023年6月頃、一部界隈で日本の住所がヤバいという話題が盛り上がりました。 詳しくはリンク先の記事を読んでいただきたいのですが、ざっくりまとめると、日本の住所には表記の揺れやローカルルールなどの問題があり、全国の住所を正規化することが困難だという内容です。 つまり収集した住所データをコンピュータの自動処理にかける前の手間がかかることを意味するのですが、このことは最近になって判明したことではなく以前からよく知られたことでした。 しかし2023年6月のTwitterでの発言をきっかけに盛り上がりました。 というわけで、ジオコーディングをサービスとして行う民間企業も以前からあったのですが、住所正規化に商機を見出したのか、この頃から各社とも積極的にサービスをPRしているような気もしています。
ただしこれらのサービスの多くは、主に法人を対象としたサービスとなっており、個人として利用することのハードルは、費用も含めて、かなり高いです。
外部のサービスとのデータのやり取りに、API1を使う方法と、ファイルを使う方法があります。
Rから外部のAPIを利用するパッケージを利用することで、簡単にジオコーディングを行うことができます。
また、外部サービスによるジオコードの結果をRに読み込み， 空間データとして分析に使用する方法も併せて紹介します。
APIを利用する方法では、tidygeocoderパッケージを利用する方法を解説します。 ファイルを利用する方法では、「地理院地図」、「CSVアドレスマッチングサービス（東京大学）」、「Google Earth Pro」を使ったジオコーディングの方法を紹介します。
使用するデータ
今回の演習では、佐賀県内にある「スターバックスコーヒー」の店舗の住所を使います。 はじめに、以下のCSVファイルをダウンロード（右クリックから[名前をつけてリンクを保存…]を選択）して、ダウンロードしたファイルをプロジェクトのdataフォルダに移動してください。
これは、佐賀県内のスターバックス店舗の住所一覧を、公式ウェブサイトの情報をもとに作成したものです。 ファイルの中身は、11店舗の店名と住所だけです。
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パッケージをロードしましょう。
library(tidygeocoder)使い方は極めて簡単です。geo関数に住所の文字列を渡すだけです。 method引数で使用するAPIを指定するのですが、ここではarcgisを使っています。
geo('佐賀県佐賀市本庄町1', method = 'arcgis')Passing 1 address to the ArcGIS single address geocoderQuery completed in: 1.3 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('福岡県北九州市小倉北区大手町11-4',
'佐賀県佐賀市本庄町1')
geo_addr <- geo(address, method = 'arcgis')Passing 2 addresses to the ArcGIS single address geocoderQuery completed in: 0.6 secondsgeo_addr# A tibble: 2 × 3
address lat long
<chr> <dbl> <dbl>
1 福岡県北九州市小倉北区大手町11-4 33.9 131.
2 佐賀県佐賀市本庄町1 33.2 130.ベクトルの長さを行数に持つデータフレームを返します。
ただしこのままでは、ただの経度と緯度の数値データからなるデータフレームなので、これを空間データに変換する必要があります。 先ほど100個のランダムポイントデータでやったのと同じ処理が必要です。
geo_addr <- geo_addr |>
st_as_sf(coords = c("long", "lat"), crs = 4326)sf::st_as_sf関数でデータフレームをsfのポイントデータに変換します。 coords引数で、データフレームのどの列を座標データに利用するかを指定し、 CRSはWGS84（EPSGコード：4326）を指定しています。
これを地図に表示してみましょう。
mapview(geo_addr) とても簡単ですね。
データフレームのバッチ処理によるジオコーディング
住所の文字列や、そのベクトルを用意すれば、ジオコードが簡単にできることがわかりました。 tidygeocoderパッケージを使えば、データフレームに対して一括してジオコードを行うこともできます。 先ほどダウンロードしたStarbucks_saga.csvをread_csv関数を使って読み込みましょう。CSVファイルのファイルエンコーディングを、RのデフォルトであるUTF-8ではなく、Shift-JISにしてあるので（Excelで開きやすいように）、locale関数を使ってShift-JISであることをread_csv関数に教えています。 店舗名（列名：名前）と住所の2列からなるデータフレームで、佐賀県内の11店舗のデータが入っていることがわかります。
starbucks <-
read_csv('data/Starbucks_saga.csv',
locale = locale(encoding = 'sjis'))Rows: 11 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: 11 × 2
名前 住所
<chr> <chr>
1 佐賀南バイパス店 佐賀県佐賀市本庄町袋306-6
2 佐賀大学通り店 佐賀県佐賀市与賀町70-1
3 鳥栖プレミアム・アウトレット店 佐賀県鳥栖市弥生が丘8-1
4 TSUTAYA鳥栖店 佐賀県鳥栖市本鳥栖町537-1
5 鳥栖蔵上町店 佐賀県鳥栖市蔵上町662-5
6 基山パーキングエリア（上り線）店 佐賀県三養基郡基山町小倉2097-1
7 蔦屋書店武雄市図書館店 佐賀県武雄市武雄町武雄5304-1
8 佐賀武雄店 佐賀県武雄市武雄町昭和277
9 唐津店 佐賀県唐津市和多田西山1-46-1
10 金立サービスエリア（下り線）店 佐賀県佐賀市金立町金立1197-57
11 ゆめタウン佐賀店 佐賀県佐賀市兵庫北5-14-1 データフレームに対してジオコーディングを行うにはgeocode関数を使います。
starbucks_tg <- starbucks |> geocode(住所, method = 'arcgis')Passing 11 addresses to the ArcGIS single address geocoderQuery completed in: 6.7 secondsstarbucks_tg# A tibble: 11 × 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 蔦屋書店武雄市図書館店 佐賀県武雄市武雄町武雄5304-1 33.2 130.
8 佐賀武雄店 佐賀県武雄市武雄町昭和277 33.2 130.
9 唐津店 佐賀県唐津市和多田西山1-46-1 33.4 130.
10 金立サービスエリア（下り線）店 佐賀県佐賀市金立町金立1197-57 33.3 130.
11 ゆめタウン佐賀店 佐賀県佐賀市兵庫北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)
mapview(starbucks_tg)地理院地図によるジオコーディング
地理院地図とは
地理院地図は、国土地理院が提供する、地図や空中写真などが閲覧できるWebサービスです。 公式サイトには、以下のように説明されています。
地理院地図とは、地形図、写真、標高、地形分類、災害情報など、国土地理院が捉えた日本の国土の様子を発信するウェブ地図です。3Dで見ることもできます。また、地形断面図の作成や新旧の写真を比較する機能なども備えています。
地理院地図によるジオコーディング
地理院地図をブラウザで開いてください。 次のような画面が表示されると思います。
このブラウザのウィンドウに、CSVファイル（Starbucks_saga.csv）をドラッグ＆ドロップします（ブラウザにファイルをドラッグ＆ドロップするという経験はあまりないと思いますので、最初は混乱するかもしれませんが、大丈夫です）。 すると、「CSVファイル読込」というウィンドウが開きますので、プルダウンメニューから[住所]を選択し、 [上記の内容で読込開始]ボタンをクリックします。
すると、佐賀県内のスターバックスの位置が、赤い丸のアイコン🔴で表示されたのがわかると思います。
住所データが、地図上に表示されたということは、その場所の緯度経度のデータがわかった、つまりジオコードされたことを意味します。
それでは、ジオコードされた緯度経度データをRに取り込む方法について説明します。 地理院地図では、ジオコードされたデータを、KMLファイルまたはGeoJSONファイルとして保存することが可能です。 ここでは、GeoJSONファイルとして保存する方法を説明します。
画面右上の[ツール]をクリックすると、画面右側にメニューが表示されますので、[作図・ファイル]をクリックします。
「作図・ファイル」というウィンドウが開きますので、[Starbucks_saga.csv]にチェックがついていることを確認したら、保存ボタン💾をクリックします。
ファイル形式を選択するウィンドウが開きますので、[GeoJSON形式]の左のラジオボタンを選択し、 下の[上記の内容で保存]をクリックします。
[保存]ボタンをクリックします。
保存したGeoJSONファイルを「Starbucks_saga.geojson」という名前に変更して、 プロジェクトのdataフォルダに移動しておきましょう。
結果の表示
重要
地理院地図によるジオコーディングの演習を行わなかった場合は、以下のリンクからGeoJSONファイルをダウンロードし、dataフォルダに移動してください。
GeoJSONファイルを、st_read関数で読み込みます。
starbucks_gsi <- st_read("data/Starbucks_saga.geojson")Reading layer `Starbucks_saga' from data source
`/Users/kzktmr/Documents/LectureNotes/data/Starbucks_saga.geojson'
using driver `GeoJSON'
Simple feature collection with 11 features and 6 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 129.9818 ymin: 33.1889 xmax: 130.5298 ymax: 33.44066
Geodetic CRS: WGS 84地理院地図で保存したGeoJSONには、地理院地図で使うためのデータも入っています。
glimpse(starbucks_gsi)Rows: 11
Columns: 7
$ name <chr> "佐賀南バイパス店", "佐賀大学通り店", "鳥栖プレミアム・…
$ 住所 <chr> "佐賀県佐賀市本庄町袋306-6", "佐賀県佐賀市与賀町70-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)mapviewを使って、地図に表示してみましょう。
mapview(starbucks_gsi)CSVアドレスマッチングサービスによるジオコーディング
CSVアドレスマッチングサービスとは
東京大学空間情報科学研究センターが提供するWebサービスで、住所情報を含むCSVファイルをアップロードすると、そのファイルに緯度経度の情報を付加してダウンロードできます。
公式サイトによる説明は、以下の通りです。
本サービスは、住所・地名フィールドを含むCSV形式データにアドレスマッチング処理を行い、緯度経度または公共測量座標系の座標値を追加します。本サービスで座標値を付加したファイルをGISで読み込めば、下のような地図を作成することもできますし、様々な空間解析を行なうこともできます。
CSVアドレスマッチングサービスによるジオコーディング
CSVアドレスマッチングサービスのサイトをブラウザで開いてください。
このページに、CSVアドレスマッチングサービスの詳しい使い方が書いてありますので、 時間のある時に目を通してみてください。
今日は、[こちらから]をクリックして、すぐにサービスを利用することにします。
対象範囲を選択します。通常は[全国街区レベル（経緯度・世界測地系）]を選んでおけばよいですが、 今回は地点が全て佐賀県内にあることがわかっていますので、[佐賀県街区レベル（経緯度・世界測地系）]を選んでみます（「経緯度」と「公共測量座標系」の違い、「世界測地系」と「旧測地系」の違いはわかりますよね？）。
住所を含むカラム番号は、今回使用するStarbucks_saga.csvでは、2列目に住所のデータがあるので、[2]と入力します。他の項目はそのままでよいです。
変換したいファイル名の[ファイルを選択]ボタンをクリックします。 ファイル選択画面が開きますので、dataフォルダにある[Starbucks_saga.csv]を選択して、[アップロード]ボタンをクリックします。
変換したいファイル名に[Starbucks_saga.csv]と書かれているのを確認したら、[送信]ボタンをクリックします。
すると、ファイルがアップロードされ、すぐにCSVファイルのダウンロードが始まると思います。 CSVファイルがダウンロードされたら、ファイル名を[Starbucks_saga_geocoded.csv]など別の名前に変えて、プロジェクトのdataフォルダに移動しましょう。
このCSVファイルをExcelで開いてみると、アップロードしたファイルに比べ、いくつかの列が追加されていることがわかると思います。 このうち、fXが経度、fYが緯度の情報を表しています。
結果の表示
重要
CSVアドレスマッチングサービスによるジオコーディングの演習を行わなかった場合は、以下のリンクからCSVファイルをダウンロードし、dataフォルダに移動してください。
この経度・緯度の情報を使って、Rでポイントデータを作成しましょう。
まずは、このCSVファイルをread_csv関数を使って読み込みます。 CSVファイルの文字コードが「sjis」なので、文字化けを防ぐためにlocale引数にその情報をlocale関数を使って与えています。
starbucks_csv <-
read_csv("data/Starbucks_saga_geocoded.csv",
locale = locale(encoding = "sjis"))Rows: 11 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (3): 名前, 住所, LocName
dbl (4): fX, fY, iConf, iLvl
ℹ 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_ut <- starbucks_csv |>
st_as_sf(coords = c("fX", "fY"), crs = 6668)
glimpse(starbucks_ut)Rows: 11
Columns: 6
$ 名前 <chr> "佐賀南バイパス店", "佐賀大学通り店", "鳥栖プレミアム・アウト…
$ 住所 <chr> "佐賀県佐賀市本庄町袋306-6", "佐賀県佐賀市与賀町70-1", "佐賀…
$ LocName <chr> "佐賀県/佐賀市/本庄町/袋/３０６番地", "佐賀県/佐賀市/与賀町",…
$ iConf <dbl> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5
$ iLvl <dbl> 7, 5, 7, 7, 7, 7, 7, 7, 7, 7, 7
$ geometry <POINT [°]> POINT (130.299 33.23963), POINT (130.2923 33.24775), POINT (1…mapviewを使って、地図に表示してみましょう。
mapview(starbucks_ut)Google Earth Proによるジオコーディング
Google Earth Proとは
Google Earthは、Googleがインターネットを前提として開発した「バーチャル地球儀システム」です。 世界中の衛星写真を、まるで地球儀を回しているかのように閲覧することができます。 Google Earthには、Google Earth（ウェブ版）とGoogle Earth Pro（デスクトップ版）があり、Google Earth Proのほうが高機能です。
Google Earth Proの入手とインストール
Google Earth Proは公式ウェブサイトから入手することができます。
[Earthのバージョン]をクリックしてください。
[Google Earthプロ（パソコン用）]をクリックしてください。
[Earthプロ（パソコン用）をダウンロード]をクリックしてください。
[同意してダウンロード]をクリックしてください。
ダウンロードが完了したら、インストールしてください。
Google Earth Proによるジオコーディング
Google Earth Proを起動してください。 下のような画面が表示されると思います。
CSVファイルを読み込むために、[ファイル]→[インポート]をクリックします。
ファイル選択ウィンドウが開きますので、先ほど[Starbucks_saga.csv]ファイルを保存したフォルダに移動します （ファイル候補が表示されない場合、右下のファイル形式を[Generic Text (* txt * csv)]に変更してください）。
すると[Starbucks_saga.csv]が表示されますので、それを選択して、右下の[開く]をクリックします。
すると、[データのインポートウィザード]が開きます。 読み込むファイルはCSV（カンマ区切り値）なので、 [フィールドタイプ]が[区切り]、[区切り]が[コンマ]になっていることを確認したら、右下の[次へ]をクリックします。
[このデータセットには練度/経度情報は含まれていませんが、番地は含まれています]にチェックを入れて、右下の[次へ]をクリックします。
[このデータセットには住所フィールドが1つあります。]のラジオボタンを選択し、 [住所フィールド]のプルダウメニューで、住所が入っている列名である[住所]を選択します。 そして、右下の[次へ]をクリックします。
フィールドタイプの指定（オプション）は、そのままで構わないので、右下の[完了]をクリックします。
ジオコードがうまくいかない住所があった場合、次のようなウィンドウが表示されることがあります。 この場合、[もしかして]をクリックし、ジオコード結果を修復することができる場合があります。
今回の場合は、見事「スターバックスコーヒー佐賀南バイパス店」を提案してくれましたので、プルダウンメニューからこれを選択し、[OK]ボタンをクリックします。
修復に成功していることを確認し、[OK]ボタンをクリックします。
取得したアイテムにスタイルデンプレートを適用しますか？と聞いてきますので、「いいえ(N)]をクリックします。
画面左の[場所]パネルにある[Starbucks_saga.csv]のチェックボックスにチェックを入れると、 スターバックス店舗の位置と店名が地図に表示されるはずです。
住所データが、地図上に表示されたということは、その場所の緯度経度のデータがわかった、つまりジオコードされたことを意味します。
それでは、ジオコードされた緯度経度データをRに取り込む方法について説明します。 Google Earth Proでは、ジオコードされたデータを、KMLファイルとして保存することが可能です。
先ほどチェックを入れた[Starbucks_saga.csv]を右クリックし、出てきたメニューから[名前を付けて場所を保存]をクリックします。
KMLファイルをdataフォルダに保存します。 プロジェクトのdataフォルダに移動し、[ファイルの種類(T):]のプルダウンメニューから[Kml(*.kml)]を選択し、[保存(S)]ボタンをクリックすれば、KMLファイルを保存することができます。
結果の表示
ノート
Google Earthの演習を行わなかった場合は、以下のリンクからKMLファイルをダウンロードし、dataフォルダに移動してください。
KMLファイルを読み込みます。st_read関数でOKです。
starbucks_ge <- st_read("data/Starbucks_saga.kml")Reading layer `Starbucks_saga' from data source
`/Users/kzktmr/Documents/LectureNotes/data/Starbucks_saga.kml'
using driver `KML'
Simple feature collection with 11 features and 2 fields
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: 129.982 ymin: 33.1888 xmax: 130.5296 ymax: 33.4408
z_range: zmin: 0 zmax: 0
Geodetic CRS: WGS 84これを見ると、DimensionがXYZとなっており、3次元のポイントデータになっていることがわかります。 しかし実際には、z_rangeを見ると、zminもzmaxも0で、高さの情報は持っていないので、st_zm関数を使って、Z次元を落としてしまいます。
starbucks_ge <- st_zm(starbucks_ge)
glimpse(starbucks_ge)Rows: 11
Columns: 3
$ Name <chr> "佐賀南バイパス店", "佐賀大学通り店", "鳥栖プレミアム・ア…
$ Description <chr> "", "", "", "", "", "", "", "", "", "", ""
$ geometry <POINT [°]> POINT (130.299 33.23874), POINT (130.2921 33.24846), POINT…mapviewを使って、地図に表示してみましょう。
mapview(starbucks_ge)ジオコーディング結果の比較
今回の演習で使用した4つのジオコーディング方法の結果を比較してみましょう。 3つの方法で作ったデータを地図に重ねて表示してみます。
mapview(starbucks_tg, col.regions = "red") +
mapview(starbucks_gsi, col.regions = "green") +
mapview(starbucks_ut, col.regions = "orange") +
mapview(starbucks_ge, col.regions = "blue")これをみると、「金立サービスエリア（下り線）店」の場所を除けば、4つのデータに大きな違いはないことがわかると思います（金立サービスエリア（下り線）店の位置については、Google Earth Proによる結果が正確のようです）。 また、地理院地図によるジオコーディングの結果と、CSVアドレスマッチングサービスを利用したジオコーディングの結果はほとんど一致しています。 地理院地図のジオコーディングは、内部的に東京大学空間情報科学研究センターが提供するWebサービスを利用しているため、両者の結果が一致するのだと考えられます。
郵便番号のジオコーディング
アンケート調査などで、居住地や勤務先の住所を尋ねるのではなく、郵便番号を尋ねることがあります。これは、有効回答率を高めるために、住所よりも回答しやすいと思われる郵便番号を尋ねることで、回答の心理的な障壁を下げるためです。
さて、郵便番号のデータがあるときに、そのジオコーディングはどのようにすればよいでしょうか。
方法の1つ目は、tidygeocoderを使うことです。 method引数にosm（Nominatim：OpenStreetMapのデータを用いたオープンソースのジオコーディングツール）を設定すれば、郵便番号によるジオコーディングが可能です。
postalcode <- c('8030814', '8408502')
postalcode_geo <- geo(postal = postalcode, method = 'osm') |>
st_as_sf(coords = c('long', 'lat'), crs = 4326)Passing 2 addresses to the Nominatim single address geocoderQuery completed in: 2.1 secondsmapview(postalcode_geo)それ以外の方法としては、日本郵便が公開している郵便番号データを利用することが考えられます。 郵便番号データダウンロードから住所の郵便番号データをダウンロードし、郵便番号をキーにして、元データと郵便番号データを結合します。 結合したデータの住所データにジオコーディングを適用すれば、結果として、郵便番号の緯度経度を得ることができます。
脚注
Application Programming Interface。ソフトウェアとWebサービスとを繋ぐインターフェイスのこと。↩︎