parzerパッケージで多様な緯度経度の表記を十進数に変換する
信頼と実績のrOpenSciから新しいパッケージがCRANに登録されました。parzerパッケージです。 このパッケージは多種多様な緯度経度の表記形式を処理し、十進数での表記(DEG: Degree)に変換してくれるものです。 （例えば「139°44’28.8869」を「139.7414」にする）
CRANリリース情報および基本的な関数の紹介はrOpenSciのブログ知りました。ぜひこちらもご覧ください。
また、私自身も以前に同様の処理を行う方法としてこのような記事を書いています。 ですがparzerパッケージを使うとより簡単に緯度経度座標の表記を修正できます。 日本語の表記を扱う際には課題が残っているので、その対策を後述します。
まずは基本的な使い方を見ていきましょう。
install.packages("parzer")
library(parzer)
パース関数
入力された座標の値を十進数での表記に修正するパース関数はparse_*()で整備されています。
parse_lat()、parse_lon()はそれぞれlatitude(緯度)、longitude(経度)を処理します。
日本経緯度原点の座標を例にします。
x <- "E139°44’28.8869" y <- "N35°39’29.1572" parse_lon(x) ## [1] 139.7414 parse_lat(y) ## [1] 35.6581
パースした結果をマップ上で確認します。せっかくなので国土地理院の地理院タイルを背景に。
library(sf) library(leaflet) library(sfheaders)
basemap <- leaflet() %>% addTiles("http://cyberjapandata.gsi.go.jp/xyz/std/{z}/{x}/{y}.png", attribution = "<a href='http://maps.gsi.go.jp/development/ichiran.html' target='_blank'>地理院タイル</a>") sfheaders::sf_point(st_point(c(parse_lon(x), parse_lat(y)))) %>% st_set_crs(value = 6668) %>% st_transform(crs = 4326) %>% mapview::mapview(map = basemap)
しっかりと「日本経緯度原点」にポイントが落ちていますね。
parse_*()は多様な表現方法に対応しています。
coords <- c(45.23323, "40:25:6N", "40°25’5.994N") parse_lat(coords) ## [1] 45.23323 40.41833 40.41833
parse_lon_lat()は二つの引数に同じ長さの経度、緯度のベクトルを与えてデータフレーム形式で結果を返却します。
df <- data.frame( lon = x, lat = y, stringsAsFactors = FALSE) parse_lon_lat(df$lon, df$lat) ## lon lat ## 1 139 35
また度分秒の要素を分解する関数として
pz_*()があります。以下の例で度分秒それぞれの要素に分解します。
pz_degree(x) ## [1] 139 pz_minute(x) ## [1] 44 pz_second(x) ## [1] 28.90863
日本語での度分秒の処理
現在のバージョン(v0.1.0)ではUnicodeを扱う際には課題があります。
https://github.com/ropensci/parzer/issues/10 github.com
具体的には次のように日本語での「東経」「北緯」、「度分秒」を扱う際に発生する問題です。
x <- "東経139度44分28秒8869" y <- "北緯35度39分29秒1572" x %>% parse_lon() ## [1] NaN y %>% parse_lat() ## [1] NaN
うまくパースされません。これに対して簡単な方法ですが、「度」「分」「秒」を変換しておきます。そうすると正常にパースできます。
@Atsushi776 さんに効率的な置換の方法を教えてもらったので書き直しています👇
str_replace_allもいいですぞ！
— atusy (@Atsushi776) 2020年3月21日
> 'ももくり3年' %>% str_replace_all(c('もも' = '桃', 'くり' = '栗'))
[1] "桃栗3年"
library(stringr) x_res <- x %>% str_replace("東経", "E") %>% str_replace_all(c("度" = "\u00b0", "分" = "\u2019", "秒" = ".")) x_res ## [1] "E139°44’28.8869" parse_lon(x_res) ## [1] 139.7414 y %>% str_replace("北緯", "N") %>% str_replace_all(c("度" = "\u00b0", "分" = "\u2019", "秒" = ".")) %>% parse_lat() ## [1] 35.6581
さらに関数化するなら以下のような感じで。
replace_dohunbyo_kanji <- function(x) { str_replace_all(x, c("東経" = "E", "西経" = "W", "北緯" = "N", "南緯" = "S"), c("E", "W", "N", "S")) %>% str_replace_all(c("度" = "\u00b0", "分" = "\u2019", "秒" = ".")) } replace_dohunbyo_kanji(x) ## [1] "E139°44’28.8869" replace_dohunbyo_kanji(y) %>% parse_lat() ## [1] 35.6581
Enjoy!
データフレーム上の緯度と経度を空間オブジェクトに変換する-🦉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 までお願いします。
