6: ジオコーディング

更新日

2025年5月21日

空間結合

GISの重要な特徴は、複数の空間データを重ねて（オーバーレイと言ったりします）分析できることです。

ここでは、ポイントデータとポリゴンデータの空間結合の例を示します。

ポイントデータにポリゴンデータを空間結合することで、それぞれのポイントデータに対して、そのポイントデータがどのポリゴン領域に含まれるかという情報を付与することができます。例えば、ある緯度経度が与えられたときに、その緯度経度がどの行政区域に含まれているか調べることができます。

データの準備

はじめに、今日の演習で使用するパッケージをロードします。

library(tidyverse)
library(sf)
library(mapview)

ここでは、佐賀県の市町ポリゴンデータと、乱数によって発生させた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/EconomicGeography/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の散布図を描いてみましょう。

ggplot(aes(x = x, y = y), data = random_points) +
  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_points

Simple feature collection with 100 features and 0 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 129.7472 ymin: 32.96266 xmax: 130.5364 ymax: 33.61195
Geodetic CRS:  WGS 84
# A tibble: 100 × 1
              geometry
 *         <POINT [°]>
 1  (130.2324 33.0636)
 2 (130.2913 33.41456)
 3 (130.3411 33.46643)
 4  (129.9261 33.5467)
 5 (130.2304 32.96341)
 6 (130.4603 33.21983)
 7 (129.8179 33.19824)
 8 (129.9092 33.15935)
 9 (130.1567 33.04395)
10 (130.4679 33.45309)
# ℹ 90 more rows

これで、100個の点からなるポイントデータrandom_pointsが作成できました。佐賀県の地図と合わせて表示してみましょう。

saga <- saga |> st_transform(6670)
random_points <- random_points |> st_transform(6670)
ggplot() +
  geom_sf(data = saga) +
  geom_sf(data = random_points)

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 [m]> POINT (-71677.42 7314.616), POINT (-65910.16 46198.09),…
$ KEN_NAME <chr> NA, "佐賀県", NA, NA, NA, NA, "佐賀県", NA, "佐賀県", NA, NA, NA, "佐賀県…
$ GST_NAME <chr> NA, "佐賀市", NA, NA, NA, NA, "西松浦郡", NA, "鹿島市", NA, NA, NA, "神埼…
$ CSS_NAME <chr> NA, NA, NA, NA, NA, NA, "有田町", NA, NA, NA, NA, NA, "吉野ヶ里町", N…
$ KEY_CODE <chr> NA, "412015", NA, NA, NA, NA, "414018", NA, "412074", NA, NA,…

これをみると、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: 47
Columns: 5
$ geometry <POINT [m]> POINT (-65910.16 46198.09), POINT (-110211.3 22606.02),…
$ KEN_NAME <chr> "佐賀県", "佐賀県", "佐賀県", "佐賀県", "佐賀県", "佐賀県", "佐賀県", "佐賀県", "佐賀県"…
$ GST_NAME <chr> "佐賀市", "西松浦郡", "鹿島市", "神埼郡", "唐津市", "藤津郡", "唐津市", "武雄市", "武雄市…
$ CSS_NAME <chr> NA, "有田町", NA, "吉野ヶ里町", NA, "太良町", NA, NA, NA, NA, NA, NA, NA…
$ KEY_CODE <chr> "412015", "414018", "412074", "413275", "412023", "414417", "…

地図に表示して、佐賀県に含まれない点が取り除かれているかどうか確認しましょう。

ggplot() +
  geom_sf(data = saga) +
  geom_sf(data = random_points_joined)

乱数から作ったポイントデータに、ポリゴンデータの持っていた情報がただしく結合されていることを確認しましょう。 filter関数を使って、佐賀市に含まれるポイントだけを表示してみると、視覚的に確認することができます。

ggplot() +
  geom_sf(data = saga) +
  geom_sf(data = filter(random_points_joined, GST_NAME == "佐賀市"))

今ここで行った操作は、ポイントデータに市町のデータを結合するものだったのですが、これは言い換えると、ポイントデータに「佐賀県佐賀市」といった住所を紐づけたことに他なりません。つまり、緯度経度のデータに、住所という属性データを結合したというわけです。

今回は市町のポリゴンデータを使いましたが、これをもっと細かいポリゴンデータ、例えば小地域データ（何丁目あるいは大字）や街区レベル（何番何号）のポリゴンデータを使えば、緯度経度を住所に変換できることがわかると思います。

ジオコーディング

ジオコーディングとは

繰り返しになりますが、先ほどの空間結合の例は、緯度経度の情報を住所に変換する基本的な考え方の説明になっています。

それでは、その逆に、住所を緯度経度に変換することは可能でしょうか？答えは、住所と緯度経度を含む空間データがあれば、可能だといえます。実際に地域分析を行う場面においては、位置データとして住所情報のみが得られる場合が多く、これを分析する際に、住所を緯度経度に変換する必要に迫られる事態に頻繁に遭遇します。

この、住所を緯度経度に変換することを、ジオコーディング（Geocoding）と呼びます。そして緯度経度を住所に変換することは、逆ジオコーディング（Reverse Geocoding）と呼びます。

Rのみでのジオコーディングはかなり煩雑なので，ここでは，外部のサーピスやアプリを使ったいくつかのジオコーディングの方法を紹介します。

ノート

2023年6月頃、一部界隈で日本の住所がヤバいという話題が盛り上がりました。詳しくはリンク先の記事を読んでいただきたいのですが、ざっくりまとめると、日本の住所には表記の揺れやローカルルールなどの問題があり、全国の住所を正規化することが困難だという内容です。つまり収集した住所データをコンピュータの自動処理にかける前の手間がかかることを意味するのですが、このことは最近になって判明したことではなく以前からよく知られたことでした。しかし2023年6月のTwitterでの発言をきっかけに盛り上がりました。というわけで、ジオコーディングをサービスとして行う民間企業も以前からあったのですが、住所正規化に商機を見出したのか、この頃から各社とも積極的にサービスをPRしているような気もしています。

ただしこれらのサービスの多くは、主に法人を対象としたサービスとなっており、個人として利用することのハードルは、費用も含めて、かなり高いです。

外部のサービスとのデータのやり取りに、API¹を使う方法と、ファイルを使う方法があります。

Rから外部のAPIを利用するパッケージを利用することで、簡単にジオコーディングを行うことができます。

また、外部サービスによるジオコードの結果をRに読み込み，空間データとして分析に使用する方法も併せて紹介します。

APIを利用する方法では、tidygeocoderパッケージを利用する方法を解説します。ファイルを利用する方法では、国土地理院のWeb地図「地理院地図」を使ったジオコーディングの方法を紹介します。

使用するデータ

今回の演習では、佐賀県内にある「スターバックスコーヒー」の店舗の住所を使います。はじめに、以下のCSVファイルをダウンロード（右クリックから[名前をつけてリンクを保存…]を選択）して、ダウンロードしたファイルをプロジェクトのdataフォルダに移動してください。

Starbucks_saga.csv

これは、佐賀県内のスターバックス店舗の住所一覧を、公式ウェブサイトの情報をもとに作成したものです。ファイルの中身は、12店舗の店名と住所だけです。

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 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: 2.1 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関数は、入力した住所のベクトルの長さと同じ行数を持つデータフレームを返します。

ただしこのままでは、ただの経度と緯度の数値データからなるデータフレームなので、これを空間データに変換する必要があります。先ほど100個のランダムポイントデータでやったのと同じ処理が必要です。

geo_addr <- 
  geo_addr |> 
  st_as_sf(coords = c("long", "lat"), crs = 4326) |> 
  st_transform(6670)
geo_addr

Simple feature collection with 3 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -103776.7 ymin: 20186.79 xmax: -66045.32 ymax: 31770.23
Projected CRS: JGD2011 / Japan Plane Rectangular CS II
# A tibble: 3 × 2
  address                                      geometry
* <chr>                                     <POINT [m]>
1 佐賀県佐賀市本庄町1              (-66045.32 27081.17)
2 佐賀県佐賀市鍋島5丁目1-1         (-68262.33 31770.23)
3 佐賀県西松浦郡有田町大野乙2441-1 (-103776.7 20186.79)

これを地図に表示してみましょう。

ggplot() +
  geom_sf(data = saga) +
  geom_sf(data = geo_addr)

とても簡単ですね。

データフレームのバッチ処理によるジオコーディング

住所の文字列や、そのベクトルを用意すれば、ジオコードが簡単にできることがわかりました。 tidygeocoderパッケージを使えば、データフレームに対して一括してジオコードを行うこともできます。先ほどダウンロードしたStarbucks_saga.csvをread_csv関数を使って読み込みましょう。CSVファイルのファイルエンコーディングを、RのデフォルトであるUTF-8ではなく、Shift-JISにしてあるので（Excelで開きやすいように）、locale関数を使ってShift-JISであることをread_csv関数に教えています。

starbucks <- 
  read_csv('data/Starbucks_saga.csv',
           locale = locale(encoding = 'sjis'))

Rows: 12 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: 12 × 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 佐賀兵庫店                       佐賀県佐賀市兵庫町藤木1192-1  
12 ゆめタウン佐賀店                 佐賀県佐賀市兵庫北5-14-1

店舗名（列名：名前）と住所の2列からなるデータフレームで、佐賀県内のR nrow(starbucks)店舗のデータが入っていることがわかります。データフレームに対してジオコーディングを行うにはgeocode関数を使います。

starbucks_tg <- starbucks |> geocode(住所, method = 'arcgis')

Passing 12 addresses to the ArcGIS single address geocoder

Query completed in: 7.2 seconds

starbucks_tg

# A tibble: 12 × 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 佐賀兵庫店                       佐賀県佐賀市兵庫町藤木1192-1    33.3  130.
12 ゆめタウン佐賀店                 佐賀県佐賀市兵庫北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) |> 
  st_transform(6670)
ggplot() +
  geom_sf(data = saga) +
  geom_sf(data = starbucks_tg)

地理院地図によるジオコーディング

地理院地図とは

地理院地図は、国土地理院が提供する、地図や空中写真などが閲覧できるWebサービスです。公式サイトには、以下のように説明されています。

地理院地図とは、地形図、写真、標高、地形分類、災害情報など、国土地理院が捉えた日本の国土の様子を発信するウェブ地図です。3Dで見ることもできます。また、地形断面図の作成や新旧の写真を比較する機能なども備えています。

地理院地図によるジオコーディング

地理院地図をブラウザで開いてください。次のような画面が表示されると思います。

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

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

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

住所データが、地図上に表示されたということは、その場所の緯度経度のデータがわかった、つまりジオコードされたことを意味します。

それでは、ジオコードされた緯度経度データをRに取り込む方法について説明します。地理院地図では、ジオコードされたデータを、KMLファイルまたはGeoJSONファイルとして保存することが可能です。ここでは、GeoJSONファイルとして保存する方法を説明します。

画面右上の[ツール]をクリックすると、画面右側にメニューが表示されますので、[作図・ファイル]をクリックします。

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

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

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

保存したGeoJSONファイルをプロジェクトのdataフォルダに移動しましょう。

結果の表示

重要

地理院地図によるジオコーディングの演習を行わなかった場合は、以下のリンクからGeoJSONファイルをダウンロードし、dataフォルダに移動してください。

Starbucks_saga.geojson

GeoJSONファイルを、st_read関数で読み込みます。

starbucks_gsi <- 
  st_read("data/Starbucks_saga.geojson") |> 
  st_transform(6670)

Reading layer `Starbucks_saga' from data source 
  `/Users/kzktmr/Documents/LectureNotes/EconomicGeography/data/Starbucks_saga.geojson' 
  using driver `GeoJSON'
Simple feature collection with 12 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: 12
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 [m]> POINT (-65331.82 26792.76), POINT (-65942.85 27697…

これはRでは使用しませんので、不要なデータを取り除きましょう。ここでは、select関数を使ってみます。

starbucks_gsi <- starbucks_gsi |> 
  select(name, address = 住所, geometry)
glimpse(starbucks_gsi)

Rows: 12
Columns: 3
$ name     <chr> "佐賀南バイパス店", "佐賀大学通り店", "鳥栖プレミアム・アウトレット店", "TSUTAYA鳥栖店", "鳥栖蔵上…
$ address  <chr> "佐賀県佐賀市本庄町袋306-6", "佐賀県佐賀市与賀町70-1", "佐賀県鳥栖市弥生が丘8-1", "佐賀県鳥栖市本…
$ geometry <POINT [m]> POINT (-65331.82 26792.76), POINT (-65942.85 27697.21),…

それでは、地図に表示してみましょう。

ggplot() + 
  geom_sf(data = saga) +
  geom_sf(data = starbucks_gsi)

ジオコーディング結果の比較

今回の演習で使用した2つのジオコーディング方法の結果を比較してみましょう。 2つの方法で作ったデータを地図に重ねて表示してみます。ここではmapviewを使いましょう。

mapview(starbucks_tg, col.regions = "red") +
  mapview(starbucks_gsi, col.regions = "green")

これをみると、「金立サービスエリア（下り線）店」の場所を除けば、2つのジオコーディング結果に大きな違いはないことがわかると思います（金立サービスエリア（下り線）店の位置については、tinygeocoderを使った結果が正確のようです）。

郵便番号のジオコーディング

アンケート調査などで、居住地や勤務先の住所を尋ねるのではなく、郵便番号を尋ねることがあります。これは、有効回答率を高めるために、住所よりも回答しやすいと思われる郵便番号を尋ねることで、回答の心理的な障壁を下げるためです。

さて、郵便番号のデータがあるときに、そのジオコーディングはどのようにすればよいでしょうか。

方法の1つ目は、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) |> 
  st_transform(6670)
mapview(postalcode_geo)

それ以外の方法としては、日本郵便が公開している郵便番号データを利用することが考えられます。郵便番号データダウンロードから住所の郵便番号データをダウンロードし、郵便番号をキーにして、元データと郵便番号データを結合します。結合したデータの住所データにジオコーディングを適用すれば、結果として、郵便番号の緯度経度を得ることができます。

saga_pc <- 
  read_csv(
    file = "data/41saga.zip", 
    col_names = FALSE, 
    locale = locale(encoding = "sjis")
  ) |> 
  select(
    postalcode = X3,
    pref = X7,
    city = X8,
    town = X9
  ) |> 
  mutate(
    postalcode = as.character(postalcode),
    address = paste0(pref, city, town)
  )

Rows: 872 Columns: 15
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (6): X4, X5, X6, X7, X8, X9
dbl (9): X1, X2, X3, X10, X11, X12, X13, X14, X15

ℹ 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.

saga_pc

# A tibble: 872 × 5
   postalcode pref   city   town                 address                        
   <chr>      <chr>  <chr>  <chr>                <chr>                          
 1 8490000    佐賀県 佐賀市 以下に掲載がない場合 佐賀県佐賀市以下に掲載がない場合……
 2 8400812    佐賀県 佐賀市 愛敬町               佐賀県佐賀市愛敬町             
 3 8400042    佐賀県 佐賀市 赤松町               佐賀県佐賀市赤松町             
 4 8400053    佐賀県 佐賀市 朝日町               佐賀県佐賀市朝日町             
 5 8400844    佐賀県 佐賀市 伊勢町               佐賀県佐賀市伊勢町             
 6 8400052    佐賀県 佐賀市 今宿町               佐賀県佐賀市今宿町             
 7 8400801    佐賀県 佐賀市 駅前中央             佐賀県佐賀市駅前中央           
 8 8400816    佐賀県 佐賀市 駅南本町             佐賀県佐賀市駅南本町           
 9 8400811    佐賀県 佐賀市 大財                 佐賀県佐賀市大財               
10 8400802    佐賀県 佐賀市 大財北町             佐賀県佐賀市大財北町           
# ℹ 862 more rows

実はread_csv関数は、zip圧縮されたファイルを読むことができます。

郵便番号データの説明ページに、郵便番号データファイルの形式が掲載されいますので、これを参考にして、3列目（郵便番号）、7列目（都道府県名）、8列目（市区町村名）、9列目（町域名）だけを残しています。都道府県名と市区町村名、町域名の文字列を連結して住所を作成しています。

post_df <- 
  tibble(postalcode = c('8400027', '8490937')) |> 
  left_join(saga_pc)

Joining with `by = join_by(postalcode)`

post_df

# A tibble: 2 × 5
  postalcode pref   city   town       address               
  <chr>      <chr>  <chr>  <chr>      <chr>                 
1 8400027    佐賀県 佐賀市 本庄町本庄 佐賀県佐賀市本庄町本庄
2 8490937    佐賀県 佐賀市 鍋島       佐賀県佐賀市鍋島

ここでは、事業用の個別郵便番号ではなく、住所の郵便番号のデータフレームを作成しています。 left_join関数で、入力した郵便番号に、住所のデータが連結されました。この住所をジオコードしてポイントデータに変換し、地図にしてみましょう。

post_df <- 
  post_df |> geocode(address, method = "arcgis") |> 
  st_as_sf(coords = c("long", "lat"), crs = 4326)

Passing 2 addresses to the ArcGIS single address geocoder

Query completed in: 0.3 seconds

mapview(post_df, zcol = "postalcode")

いかがでしょうか。手間はかかりますが、先ほどのtidygeocoderを使ったのと同じような結果が得られました。

このように、さまざまなデータを組み合わせて、工夫して空間データを作成することができれば、地域分析の幅が広がります。

脚注

Application Programming Interface。ソフトウェアとWebサービスとを繋ぐインターフェイスのこと。↩︎