データフレーム上の緯度と経度を空間オブジェクトに変換する-🦉sfと🐍geopandasの例-
地理空間データを取り扱う際は、はじめにgeojsonやshapeファイルで受け取ることが一般的かと思います。 あるいはポイントデータの場合には、緯度と経度の値が各列に記録されるcsvなどの表形式のテキストファイルを起点とすることもあります。
前者のような地理空間データであれば、 R等のアプリケーションで読み込めば自動的に地理空間データとみなしてくれますが、 後者の場合はそうではありません。 座標の値はあくまでも数値です。なのでこうしたデータを地理空間データとして扱えるようにするには変換作業が必要となります。
今回の記事では、RおよびPythonでのデータフレームに記録された緯度経度の列を変換し、空間オブジェクトとして扱えるようにする方法を紹介します。 空間オブジェクトの形式としてRではsf、Pythonではgeopandasを対象にします。地理空間データに変換しておくと、データの空間的な配置を可視化可能になったり空間演算が可能になり、分析の幅を広げることが期待できます。
対象データ
csvファイルに緯度と経度の値が格納されているデータの例として、国土交通省 位置参照情報ダウンロードサービス 提供の平成30年整備 茨城県つくば市の大字・町丁目レベルのデータを利用します。
ファイルをダウンロードしたら、まずはデータフレームとしてcsvファイルを読み込みましょう。
library(sf) library(readr)
df_isj08220 <- read_csv("08220_2018.csv", locale = locale(encoding = "cp932"), col_types = "___c_cdd__") %>% purrr::set_names(c("city", "street_lv1", "latitude", "longitude"))
データの一部を表示します。latitude, longitudeの列がそれぞれ独立しています。
| city | street_lv1 | latitude | longitude |
|---|---|---|---|
| つくば市 | 赤塚 | 36.04268 | 140.1232 |
| つくば市 | 明石 | 36.18207 | 140.0537 |
| つくば市 | 旭 | 36.11289 | 140.0725 |
| つくば市 | 安食 | 36.16775 | 140.0163 |
| つくば市 | あしび野 | 35.96060 | 140.1122 |
| つくば市 | 吾妻一丁目 | 36.08153 | 140.1129 |
sfパッケージ
sfパッケージでデータフレームの座標情報を参照してポイントデータに変換するにはst_as_sf()を使います。
第一引数に座標列を含む対象のデータフレーム、coords引数の値にX,Y座標となる列名を指定します。この時、座標参照系が判明している際はcrs引数によって明示しておきましょう。
sf_isj08220 <- df_isj08220 %>% sf::st_as_sf(coords = c("longitude", "latitude"), crs = 4326)
sf_isj08220 ## Simple feature collection with 299 features and 2 fields ## geometry type: POINT ## dimension: XY ## bbox: xmin: 140.0102 ymin: 35.95047 xmax: 140.1676 ymax: 36.22462 ## epsg (SRID): 4326 ## proj4string: +proj=longlat +datum=WGS84 +no_defs ## # A tibble: 299 x 3 ## city street_lv1 geometry ## <chr> <chr> <POINT [°]> ## 1 つくば市 赤塚 (140.1232 36.04268) ## 2 つくば市 明石 (140.0537 36.18207) ## 3 つくば市 旭 (140.0725 36.11289) ## 4 つくば市 安食 (140.0163 36.16775) ## 5 つくば市 あしび野 (140.1122 35.9606) ## 6 つくば市 吾妻一丁目 (140.1129 36.08153) ## 7 つくば市 吾妻二丁目 (140.1111 36.08443) ## 8 つくば市 吾妻三丁目 (140.1109 36.0889) ## 9 つくば市 吾妻四丁目 (140.1183 36.08722) ## 10 つくば市 天久保一丁目 (140.1087 36.0919) ## # … with 289 more rows
緯度経度の列を元にポイントデータが作成できました。地図上にマッピングしてみます。
mapview::mapview(sf_isj08220)
sfheadersパッケージ
続いてsfheadersパッケージを使う方法です。sf::st_as_sf()の実行結果と同じく、sfオブジェクトを作成します。
今回は対象がポイントデータなので、sf_point()のx,y引数にそれぞれ経度と緯度の列名を与えます。
この時、keep引数の値をTRUEにしておくとデータフレームの他の列が残ります。FALSEではポイントのgeometry列のみしか残らないため注意です。
また、座標参照系もこのタイミングでは指定できないため、sf::st_set_crs()で定義する必要が生じます。
df_isj08220 %>% sfheaders::sf_point(x = "longitude", y = "latitude", keep = TRUE) %>% st_set_crs(4326) %>% tibble::new_tibble(subclass = "sf", nrow = nrow(.))
sf_isj08220sfh ## Simple feature collection with 299 features and 2 fields ## geometry type: POINT ## dimension: XY ## bbox: xmin: 140.0102 ymin: 35.95047 xmax: 140.1676 ymax: 36.22462 ## z_range: zmin: NA zmax: NA ## m_range: mmin: NA mmax: NA ## epsg (SRID): 4326 ## proj4string: +proj=longlat +datum=WGS84 +no_defs ## # A tibble: 299 x 3 ## city street_lv1 geometry ## <chr> <chr> <POINT [°]> ## 1 つくば市 赤塚 (140.1232 36.04268) ## 2 つくば市 明石 (140.0537 36.18207) ## 3 つくば市 旭 (140.0725 36.11289) ## 4 つくば市 安食 (140.0163 36.16775) ## 5 つくば市 あしび野 (140.1122 35.9606) ## 6 つくば市 吾妻一丁目 (140.1129 36.08153) ## 7 つくば市 吾妻二丁目 (140.1111 36.08443) ## 8 つくば市 吾妻三丁目 (140.1109 36.0889) ## 9 つくば市 吾妻四丁目 (140.1183 36.08722) ## 10 つくば市 天久保一丁目 (140.1087 36.0919) ## # … with 289 more rows
おまけ: sfからdfに戻す
データフレーム(df)の緯度経度情報をもとにsfオブジェクトを作成する例を述べましたが、反対にsfからデータフレームへ、地理空間情報を取り除く方法も紹介します。
作成したsfオブジェクトに対してsf::st_drop_geometry()ないしsf::st_set_geometry(value = NULL)でgeometry情報を除外することでdfに変換されます。
sf_isj08220 %>% st_drop_geometry() sf_isj08220 %>% st_set_geometry(value = NULL)
この時、ポイントの情報も失われるので、データフレームの列として残したい場合は予め列を追加する必要が生じます。
sf_isj08220 %>% dplyr::mutate(coords = purrr::map(geometry, ~ as.data.frame(sf::st_coordinates(.x)))) %>% tidyr::unnest(cols = coords) %>% st_drop_geometry() ## # A tibble: 299 x 4 ## city street_lv1 X Y ## * <chr> <chr> <dbl> <dbl> ## 1 つくば市 赤塚 140. 36.0 ## 2 つくば市 明石 140. 36.2 ## 3 つくば市 旭 140. 36.1 ## 4 つくば市 安食 140. 36.2 ## 5 つくば市 あしび野 140. 36.0 ## 6 つくば市 吾妻一丁目 140. 36.1 ## 7 つくば市 吾妻二丁目 140. 36.1 ## 8 つくば市 吾妻三丁目 140. 36.1 ## 9 つくば市 吾妻四丁目 140. 36.1 ## 10 つくば市 天久保一丁目 140. 36.1 ## # … with 289 more rows
Python
import pandas as pd import geopandas
df = pd.read_csv("08220_2018.csv", encoding = "cp932", usecols = [3, 5, 6, 7], dtype = {'市区町村名': str, '大字町丁目名': str, '緯度': float, '経度': float}) df = df.rename(columns = {'市区町村名': 'city', '大字町丁目名': 'street_lv1', '緯度': 'latitude', '経度': 'longitude'})
geopandas
gdf = geopandas.GeoDataFrame(df, geometry = geopandas.points_from_xy(df.longitude, df.latitude))
gdf.head(10)
| city | street_lv1 | longitude | latitude | geometry |
|---|---|---|---|---|
| つくば市 | 赤塚 | 140.1 | 36.04 | POINT (140.12315 36.042677) |
| つくば市 | 明石 | 140.1 | 36.18 | POINT (140.053737 36.182065) |
| つくば市 | 旭 | 140.1 | 36.11 | POINT (140.072495 36.11289) |
| つくば市 | 安食 | 140.0 | 36.17 | POINT (140.016345 36.167751) |
| つくば市 | あしび野 | 140.1 | 35.96 | POINT (140.112185 35.960604) |
| つくば市 | 吾妻一丁目 | 140.1 | 36.08 | POINT (140.112938 36.081534) |
| つくば市 | 吾妻二丁目 | 140.1 | 36.08 | POINT (140.111144 36.084425) |
| つくば市 | 吾妻三丁目 | 140.1 | 36.09 | POINT (140.11092 36.088898) |
| つくば市 | 吾妻四丁目 | 140.1 | 36.09 | POINT (140.118261 36.087219) |
| つくば市 | 天久保一丁目 | 140.1 | 36.09 | POINT (140.108708 36.091901) |
せっかくなのでこちらもマッピング。
import geopatra m = gdf.folium.plot(zoom = 10) m
いいですね！Enjoy!
郵便番号データの操作と祝日判定の機能を追加したzipangu v0.2.0を公開
昨年バージョン0.1.0をリリースした、日本人が扱う住所や年号、漢数字のデータ操作を楽にするRパッケージ、{zipangu}に新しい機能を追加し、バージョン0.2.0としてCRANに登録しました。この記事では0.2.0で扱える新機能を紹介します。表題の通り、郵便番号データの操作と、祝日の判定機能がメインです。
https://cran.r-project.org/web/packages/zipangu/news/news.html
library(zipangu) library(dplyr, warn.conflicts = FALSE)
packageVersion("zipangu") ## [1] '0.2.0'
1. 郵便番号
日本郵便が住所別の郵便番号データを公開しています。このファイルをRに取り込む関数を用意しました。郵便番号ファイルは頻繁に更新されるため、ローカルでのパス指定だけでなくURLによるファイル読み込みも可能となっています。
公開されるファイルには、住所の表記に関していくつかの種類や事業所を含めたデータがありますが、こうしたバリエーションにも対応します。具体的には、type引数に対象ファイルを元に以下の値を与えます。
oogaki: 読み仮名データの促音・拗音を小書きで表記しないものkogaki: 読み仮名データの促音・拗音を小書きで表記するものroman: 読み仮名データがローマ字jigyosyo: 事業所の個別郵便番号
# パッケージに含まれるサンプルデータを読み込みます zipcode_file <- system.file("zipcode_dummy/13TOKYO_kogaki.CSV", package = "zipangu") df_oogaki <- zipangu::read_zipcode(path = zipcode_file, type = "oogaki")
glimpse(df_oogaki) ## Observations: 1 ## Variables: 15 ## $ jis_code <chr> "13101" ## $ old_zip_code <chr> "100" ## $ zip_code <chr> "1000001" ## $ prefecture_kana <chr> "トウキヨウト" ## $ city_kana <chr> "チヨダク" ## $ street_kana <chr> "チヨダ" ## $ prefecture <chr> "東京都" ## $ city <chr> "千代田区" ## $ street <chr> "千代田" ## $ is_street_duplicate <dbl> 0 ## $ is_banchi <dbl> 0 ## $ is_cyoumoku <dbl> 0 ## $ is_zipcode_duplicate <dbl> 0 ## $ status <dbl> 0 ## $ modify_type <dbl> 0
URLを直接指定する際は type の指定が不要です。
read_zipcode("https://www.post.japanpost.jp/zipcode/dl/jigyosyo/zip/jigyosyo.zip")
read_zipcode()
では圧縮ファイルをRに読み込みますが、dl_zipcode_file()によりファイルをダウンロードしておくことも可能です。
dl_zipcode_file(path = "https://www.post.japanpost.jp/zipcode/dl/oogaki/zip/02aomori.zip")
このほか、is_zipcode()と zipcode_spacer()
を用意しました。前者は入力が7桁の郵便番号かどうかをチェック、後者は郵便番号に使われるハイフンを追加・除去するちょっとした関数です。
is_zipcode(7000027) ## [1] TRUE
is_zipcode("700-0027") ## [1] TRUE
zipcode_spacer("3050053") ## [1] "305-0053"
# ハイフンを取り除く際は remove = TRUE を指定します zipcode_spacer("305-0053", remove = TRUE) ## [1] "3050053"
郵便番号データの公開は貴重ですが「KEN_ALL.csvの闇」と噂される記録のされ方で、実際に扱う場合は処理が面倒です。こうしたデータをtidyにする試みについては前回書いたとおりです。
この課題に関して、id:yutannihilation さんにも取り込んでいただきました。
機能要望として郵便番号の検索が提案されています。どういう形でデータを整形・表示すれば良いのかについて、議論の余地が残ると判断し、完成に至っていません。
2. 日本の祝日
祝日（国民の休日）判定のための関数を用意しました。これらは{zipangu}パッケージの開発のきっかけとなる{Nippon}で提供されていた機能でもあります。最初のリリースでは実装が追いついていませんでしたので、当初の目的を果たせました（拍手）。
入力に与えた日付が祝日かを返却する is_jholiday() と その年の祝日を調べる jholiday() があります。
Nippon::is.jholiday(as.Date("2019-12-23")) ## [1] TRUE is_jholiday(date = as.Date("2019-12-23")) ## [1] FALSE
# 今年の祝日を調べる jholiday(year = 2020, lang = "jp") ## $元日 ## [1] "2020-01-01" ## ## $成人の日 ## [1] "2020-01-13" ## ## $建国記念の日 ## [1] "2020-02-11" ## ## $天皇誕生日 ## [1] "2020-02-23" ## ## $春分の日 ## [1] "2020-03-20" ## ## $昭和の日 ## [1] "2020-04-29" ## ## $憲法記念日 ## [1] "2020-05-03" ## ## $みどりの日 ## [1] "2020-05-04" ## ## $こどもの日 ## [1] "2020-05-05" ## ## $海の日 ## [1] "2020-07-23" ## ## $スポーツの日 ## [1] "2020-07-24" ## ## $山の日 ## [1] "2020-08-10" ## ## $敬老の日 ## [1] "2020-09-21" ## ## $秋分の日 ## [1] "2020-09-22" ## ## $文化の日 ## [1] "2020-11-03" ## ## $勤労感謝の日 ## [1] "2020-11-23"
# 天皇誕生日は時代に合わせて変更されます jholiday_spec(year = 1988, name = "天皇誕生日", lang = "jp") ## [1] "1988-04-29"
jholiday_spec(2018, "天皇誕生日", lang = "jp") ## [1] "2018-12-23"
jholiday_spec(2020, "天皇誕生日", lang = "jp") ## [1] "2020-02-23"
大元の祝日の一覧は内閣府が配布するcsvファイルを参照しています。2020年1月現在に決まっているものなので、将来、祝日の変更があった際はデータを更新しなくてはいけなくなりますが、アーカイブされてしまった{Nippon}で用意されていた機能を復活させたこと、天皇の即位や東京オリンピック開催に関する2019、2020年の祝日の変更に対応できた点は良かったと思います。
また、ここでの作業の副産物として「ある月の第二月曜日」などを調べるための関数 find_date_by_wday() ができました。年
(year)、月 (month)、曜日 (wday… デフォルトでは日曜日が1)、週番号 (ordinal)
を指定して実行します。
find_date_by_wday(year = 2020, month = 1, wday = 2, ordinal = 2) ## [1] "2020-01-13"
不具合や新機能の提案はGitHub issuesまたは Twitter @u_ribo までお願いします。
郵便番号データをtidyにする挑戦
日本人が頻繁に遭遇するデータ操作を効率的に行うための{zipangu}パッケージ、想定よりも多くの人が喜んでくれたようで、私としても嬉しく思っています。
はてなブログに投稿しました #はてなブログ
— Uryu Shinya (@u_ribo) 2019年12月2日
住所や年号、漢数字のデータ操作を楽にするRパッケージをCRANに登録しました - cucumber fleshhttps://t.co/5I3rntfrio
記事の最後にプロジェクトの協力者を募集したら数名からの反応があり、また新機能の要望も挙げられました。 ありがとうございます。
さて、次のリリースでは郵便番号の処理を効率的に行う機能を実装する計画でいます。 具体的には日本郵便が提供する郵便番号csvファイル（以下、郵便番号データ）をRで読み込む関数と、郵便番号の検索および住所情報を返却する機能です。
Add japan zip-code utility functions by uribo · Pull Request #6 · uribo/zipangu · GitHub
ファイルの読み込みに関してはすでにmasterブランチへマージされています。 そして郵便番号検索の方もここで読み込んだファイルを利用すれば良かろうと思っていたのですが、こんなご意見をいただきました。
まじか、Pythonに続いてRもken_ALL.csvの闇に突貫するのか… https://t.co/MKV7wjdJpa
— Aki Ariga (@chezou) 2019年12月10日
どうもこの郵便番号データには問題があるそうです。探ってみましょう。
install.packages("remotes") remotes::install_github("uribo/zipan")
library(dplyr) library(zipangu) library(stringr)
# read_zipcode() が郵便番号データを読み込むための関数です。 # 提供されている3種類(読み仮名データの促音・拗音を小書きで表記しないもの、読み仮名データの促音・拗音を小書きで表記するもの、ローマ字)の住所の表記形式、事業所のcsvファイルに対応します # path引数にzipファイルが置かれるURLまたはコンピュータ上のファイルパスを指定します df <- read_zipcode(path = "https://www.post.japanpost.jp/zipcode/dl/oogaki/zip/ken_all.zip", type = "kogaki") %>% # 市区町村コード、郵便番号、住所に関する列、重複の判定のために「丁目を有する町域の場合の表示」の列を選んでおきます select(jis_code, zip_code, prefecture, city, street, is_cyoumoku)
例えば郵便番号066-0005のレコードを検索すると次のように3件のデータが返却されます。同じ郵便番号、市町村なのになぜ？となりますが、street列で表示している町域名がおかしいことに気がつきます。
df %>% filter(zip_code == "0660005")
| jis_code | zip_code | prefecture | city | street | is_cyoumoku |
|---|---|---|---|---|---|
| 01224 | 0660005 | 北海道 | 千歳市 | 協和(88−2、271−10、343−2、404−1、427− | 0 |
| 01224 | 0660005 | 北海道 | 千歳市 | 3、431−12、443−6、608−2、641−8、814、842− | 0 |
| 01224 | 0660005 | 北海道 | 千歳市 | 5、1137−3、1392、1657、1752番地) | 0 |
なんと、これは郵便番号データ番号ファイルの仕様です。データの説明書きに次の記述があります。
全角となっている町域部分の文字数が38文字を越える場合、また半角となっているフリガナ部分の文字数が76文字を越える場合は、複数レコードに分割しています
このままでは検索用の関数を用意する際に問題になります。
また、これ以外にも数々の問題点があり、これまでに多くの方が記事にまとめてくださっています。 この行の分割をはじめとしていくつかの問題への対策方法が書かれている記事も見受けられました。
一方で、次の2点に関する具体的な処理方法については見つけることができませんでした。
藤野(400、400−2番地)など「、」で複数の住所がある –>藤野400、藤野400-2番地の行を分ける大通西(1〜19丁目)のように住所が省略される –>大通西1丁目、大通西2丁目、大通西3丁目、…大通西19丁目を独立させる
以下に示すように、元のデータはtidy1ではありません。 扱うデータがtidyであることを心がける身としては放っては置けない問題です。
df %>% filter(zip_code %in% c("0050840", "0600042"))
| jis_code | zip_code | prefecture | city | street | is_cyoumoku |
|---|---|---|---|---|---|
| 01101 | 0600042 | 北海道 | 札幌市中央区 | 大通西(1〜19丁目) | 1 |
| 01106 | 0050840 | 北海道 | 札幌市南区 | 藤野(400、400−2番地) | 0 |
このデータをtidyにするならこうかなと思います (大通西の住所は一部省略)
| jis_code | zip_code | prefecture | city | street |
|---|---|---|---|---|
| 01106 | 0050840 | 北海道 | 札幌市南区 | 藤野400 |
| 01106 | 0050840 | 北海道 | 札幌市南区 | 藤野400-2番地 |
| 01101 | 0600042 | 北海道 | 札幌市中央区 | 大通西1丁目 |
| 01101 | 0600042 | 北海道 | 札幌市中央区 | 大通西2丁目 |
そんなわけで前置きが長くなりましたが、こうした問題の解決に取り組んでいます。 いくつかの課題に関しては解決できそうと目処が立つ、一方で完璧には程遠いことを感じてきたので一旦整理しておきます。
目次
- 住所の重複と複数行に分割される問題への対処
- 「、」で区切られた住所を分割する
- 京都市内の通り名
- 複雑な町域名
- 「〜」によって省略される住所を復元する
- 残った課題
-
効果的なデータ分析を行いやすくするためのデータの持ち方を指す概念。参考: https://www.jstatsoft.org/article/view/v059i10↩
