kokudosuuchiは1年ほど前にCRANに公開したパッケージで、国土数値情報APIをRから使えます。
ここ最近、国土数値情報ウェブサイトのHTMLテーブルなんかと死闘を繰り広げていたわけですが、そろそろ「ふん、今日はこの辺で勘弁しといてやるか...（肩で息をしながら」みたいな程度には終わったので、ブログに書きます。
アップデート内容は主に以下の2つです。
getKSJData()がsf形式でデータを読み込むようになりました。- そこそこの精度でメタデータを紐づけられるようになりました。
国土数値情報ダウンロードサービスのAPIを使うための関数には変化はありません。
インストール
開発版はinstall_github()でインストールできます。
devtools::install_github("yutannihilation/kokudosuuchi")
注意点としては、次のバージョンからsfパッケージに依存するようになっています。Windowsユーザは何も考えずにインストールできますが、MacやLinuxの場合は依存ライブラリをあらかじめインストールする必要があります。 以下の指示に従ってください。
もし、sfをインストールできない環境でkokudosuuchiパッケージを使っているので困る、という方は、このブログのコメント欄とかTwitterとかでご相談ください。
利用規約
kokudosuuchiパッケージを読み込むときにも出るように、国土数値情報ダウンロードサービスは利用約款を確認の上使いましょう。
library(kokudosuuchi) #> このサービスは、「国土交通省 国土数値情報（カテゴリ名）」をもとに加工者が作成 #> 以下の国土数値情報ダウンロードサービスの利用約款をご確認の上ご利用ください： #> #> http://nlftp.mlit.go.jp/ksj/other/yakkan.html
また、データごとに利用上の注意とかがあったりもするので、各データの詳細ページについても読みましょう。
...詳細ページ？ そんなものどこにあるの？と思ったあなた、鋭い。そのことについてはあとで触れます。
getKSJData()
国土数値情報ダウンロードサービスWeb APIは、Shapefileを含むZIPファイルのURLを返してくれます。getKSJData()は、そのデータをsf形式で読み込むための関数です。
ZIPファイルのURLを探す
まずはURLを探します。台風も近づいているので、浸水想定区域を調べてみましょう。まずは「浸水想定区域」に対応する識別子（identifier）を探すためにgetKSJSummary()を使います。
getKSJSummary() %>% filter(title == "浸水想定区域") #> # A tibble: 1 x 5 #> identifier title field1 field2 areaType #> <chr> <chr> <chr> <chr> <chr> #> 1 A31 浸水想定区域 政策区域 災害・防災 3
識別子はA31のようです。これを指定してgetKSJURL()でZIPファイルのURLを取得しましょう。choose_prefecture_code()はインタラクティブに都道府県のコードを選択するための関数です。チェックボックスにチェックを入れるとその都道府県のコードを返してくれます。今回は東京都に絞り込んでみます。
urls <- getKSJURL("A31", prefCode = choose_prefecture_code()) #> Loading required package: shiny #> #> Listening on http://127.0.0.1:7648
結果は以下のようになっています。
glimpse(urls) #> Observations: 1 #> Variables: 9 #> $ identifier <chr> "A31" #> $ title <chr> "浸水想定区域" #> $ field <chr> "政策区域" #> $ year <chr> "2012" #> $ areaType <chr> "3" #> $ areaCode <chr> "13" #> $ datum <chr> "1" #> $ zipFileUrl <chr> "http://nlftp.mlit.go.jp/ksj/gml/data/A31/A31-12/A31-12_13_GML.zip" #> $ zipFileSize <chr> "2.53MB"
データを読み込むその前に
さて、これでもうZIPファイルのURLを指定すればデータを読み込めてしまうわけなんですが、このデータの利用規約は（国土数値情報全体の利用約款に加えて）データの詳細ページに書かれています。
しかし。
上の結果を見ても分かるように、なんと、データの詳細ページがどこにあるのかAPIは教えてくれません。。
そこで、getKSJData()はデータから推測して詳細ページのURLをサジェストしてくれます。
d <- getKSJData(urls$zipFileUrl, cache_dir = "cached_zip") #> #> Details about this data can be found at http://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-A31.html #>
ということで、リンク先に従って原典を以下に明示します。上に貼ったような個別の「データの使用許諾条件」はこのデータにはないようです。
河川管理者（国土交通大臣、都道府県知事）「浸水想定区域図」 ※水防法第十条第二項及び第十一条第一項に基づき指定される洪水予報河川並びに水防法第十三条に基づき指定される水位周知河川のうち、各河川管理者より資料提供を受けられたもの。
利用条件がデータによって微妙に違っていて確認が面倒なんですが、なるべく気を付けつつ使うようにしましょう。
データの読み込み
さて、さっきもう読み込んでしまいましたが、getKSJData()でデータを読み込めます。cache_dir引数は、ダウンロードしたZIPファイルを保存しておくディレクトリを指定します（指定しなければtempdir()が使われます）。
ちなみに、getKSJData()にはURLだけでなく、ダウンロードしたZIPファイルや、それを展開したディレクトリを指定することもできます。
# ZIPファイルを指定 d <- getKSJData("cached_zip/A31-12_13_GML.zip") # ディレクトリを指定 zip_exdir <- tempfile() dir.create(zip_exdir) unzip("cached_zip/A31-12_13_GML.zip", exdir = zip_exdir) d <- getKSJData(zip_exdir)
データの中身はこんな感じです。
d #> $`A31-12_13` #> Simple feature collection with 7303 features and 6 fields #> geometry type: MULTIPOLYGON #> dimension: XY #> bbox: xmin: 139.3212 ymin: 35.53492 xmax: 139.9186 ymax: 35.81771 #> epsg (SRID): NA #> proj4string: +proj=longlat +ellps=GRS80 +no_defs #> # A tibble: 7,303 x 7 #> A31_001 A31_002 A31_003 A31_004 A31_005 A31_006 geometry #> <int> <int> <chr> <chr> <chr> <chr> <S3: sfc_MULTIPOLYGON> #> 1 11 11 埼玉県 県土整備部 河川砂防課 平成19年3月27日 埼玉県告示第504号 11005 <S3: sfc_MULTIPOLYGON> #> 2 11 11 埼玉県 県土整備部 河川砂防課 平成19年3月27日 埼玉県告示第504号 11005 <S3: sfc_MULTIPOLYGON> #> 3 11 11 埼玉県 県土整備部 河川砂防課 平成19年3月27日 埼玉県告示第504号 11005 <S3: sfc_MULTIPOLYGON> #> 4 11 11 埼玉県 県土整備部 河川砂防課 平成19年3月27日 埼玉県告示第504号 11005 <S3: sfc_MULTIPOLYGON> #> 5 11 11 埼玉県 県土整備部 河川砂防課 平成19年3月27日 埼玉県告示第504号 11005 <S3: sfc_MULTIPOLYGON> #> 6 11 11 埼玉県 県土整備部 河川砂防課 平成19年3月27日 埼玉県告示第504号 11005 <S3: sfc_MULTIPOLYGON> #> 7 11 11 埼玉県 県土整備部 河川砂防課 平成19年3月27日 埼玉県告示第504号 11005 <S3: sfc_MULTIPOLYGON> #> 8 11 11 埼玉県 県土整備部 河川砂防課 平成19年3月27日 埼玉県告示第504号 11005 <S3: sfc_MULTIPOLYGON> #> 9 11 11 埼玉県 県土整備部 河川砂防課 平成19年3月27日 埼玉県告示第504号 11005 <S3: sfc_MULTIPOLYGON> #> 10 13 11 埼玉県 県土整備部 河川砂防課 平成19年3月27日 埼玉県告示第504号 11005 <S3: sfc_MULTIPOLYGON> #> # ... with 7,293 more rows #>
translateKSJData()
さて、データは読み込めましたが、上のデータのままでは列名がA31_001、A31_002といった記号だし、中身もよく分からない数字だし、ぱっと見ではよくわかりません。これらの列や数字が何を表しているかは、さきほどの詳細ページに書かれています。
例えば、A31_001の列が表す内容は以下の通りです。
「浸水深ランクコード」というリンク先には、この列の数字とそれが表す意味の対応表があります。
ということで、正攻法としてはこれをひとつひとつ紐づけていくことなんですが、それはめんどくさいので、そこそこの精度で紐づけをやってくれる関数translateKSJData()をつくりました。具体的にはこんな感じです。
d_translated <- translateKSJData(d) d_translated #> $`A31-12_13` #> Simple feature collection with 7303 features and 6 fields #> geometry type: MULTIPOLYGON #> dimension: XY #> bbox: xmin: 139.3212 ymin: 35.53492 xmax: 139.9186 ymax: 35.81771 #> epsg (SRID): NA #> proj4string: +proj=longlat +ellps=GRS80 +no_defs #> # A tibble: 7,303 x 7 #> 浸水深 作成種別 作成主体 指定年月日 告示番号 A31_006 geometry #> <chr> <chr> <chr> <chr> <chr> <chr> <S3: sfc_MULTIPOLYGON> #> 1 0～0.5ｍ未満（5段階） 埼玉県 埼玉県 県土整備部 河川砂防課 平成19年3月27日 埼玉県告示第504号 11005 <S3: sfc_MULTIPOLYGON> #> 2 0～0.5ｍ未満（5段階） 埼玉県 埼玉県 県土整備部 河川砂防課 平成19年3月27日 埼玉県告示第504号 11005 <S3: sfc_MULTIPOLYGON> #> 3 0～0.5ｍ未満（5段階） 埼玉県 埼玉県 県土整備部 河川砂防課 平成19年3月27日 埼玉県告示第504号 11005 <S3: sfc_MULTIPOLYGON> #> 4 0～0.5ｍ未満（5段階） 埼玉県 埼玉県 県土整備部 河川砂防課 平成19年3月27日 埼玉県告示第504号 11005 <S3: sfc_MULTIPOLYGON> #> 5 0～0.5ｍ未満（5段階） 埼玉県 埼玉県 県土整備部 河川砂防課 平成19年3月27日 埼玉県告示第504号 11005 <S3: sfc_MULTIPOLYGON> #> 6 0～0.5ｍ未満（5段階） 埼玉県 埼玉県 県土整備部 河川砂防課 平成19年3月27日 埼玉県告示第504号 11005 <S3: sfc_MULTIPOLYGON> #> 7 0～0.5ｍ未満（5段階） 埼玉県 埼玉県 県土整備部 河川砂防課 平成19年3月27日 埼玉県告示第504号 11005 <S3: sfc_MULTIPOLYGON> #> 8 0～0.5ｍ未満（5段階） 埼玉県 埼玉県 県土整備部 河川砂防課 平成19年3月27日 埼玉県告示第504号 11005 <S3: sfc_MULTIPOLYGON> #> 9 0～0.5ｍ未満（5段階） 埼玉県 埼玉県 県土整備部 河川砂防課 平成19年3月27日 埼玉県告示第504号 11005 <S3: sfc_MULTIPOLYGON> #> 10 1.0～2.0ｍ未満（5段階） 埼玉県 埼玉県 県土整備部 河川砂防課 平成19年3月27日 埼玉県告示第504号 11005 <S3: sfc_MULTIPOLYGON> #> # ... with 7,293 more rows #>
どうでしょうか。だいたい人間が読める名前に変換されています。
しかし、ちょっと右にスクロールしてもらうとわかりますが、A31_006という列がそのままになっています。この列は、データ詳細のページにも載っていない謎の列です。
可能な限りいい感じにメタデータを紐づけられるようにがんばったつもりなんですが、こんな感じの微妙におかしいデータがあったりして、必ず紐づけが成功するわけではありません。また、紐づけられていてもそれが正確とは限らないので（kokudosuuchiパッケージが間違っていることもあれば、元データが間違っていることもある）、紐づけた後でそれがあっているか確認をお願いします。
あと、translateKSJData()の結果がおかしいときはお知らせいただけるとうれしいです。
可視化
このデータは例えば、GitHub版のggplot2を使って以下のように可視化できます。
library(sf) d_summarised <- d_translated$`A31-12_13` %>% mutate(depth_min = readr::parse_number(stringr::str_extract(浸水深, "^[\\d\\.]+"))) %>% group_by(depth_min) %>% summarise() d_summarised #> Simple feature collection with 7 features and 1 field #> geometry type: MULTIPOLYGON #> dimension: XY #> bbox: xmin: 139.3212 ymin: 35.53492 xmax: 139.9186 ymax: 35.81771 #> epsg (SRID): NA #> proj4string: +proj=longlat +ellps=GRS80 +no_defs #> # A tibble: 7 x 2 #> depth_min geometry #> <dbl> <simple_feature> #> 1 0.0 <MULTIPOLYGON...> #> 2 0.5 <MULTIPOLYGON...> #> 3 1.0 <MULTIPOLYGON...> #> 4 2.0 <MULTIPOLYGON...> #> 5 3.0 <MULTIPOLYGON...> #> 6 4.0 <MULTIPOLYGON...> #> 7 5.0 <MULTIPOLYGON...> library(ggplot2) ggplot(d_summarised) + geom_sf(aes(fill = depth_min))
mapviewパッケージとかでもいいかもしれません（こっちはまだ使い方よく分からないんですけどこれであってます...？）
mapview::mapview(d_summarised) #> Warning message: #> In leaflet_sfc(sf::st_geometry(x), map = map, zcol = zcol, color = clrs, : #> the supplied feature layer has more points/vertices than the set threshold. #> using special rendering function, hence things may not behave as expected from a standard leaflet map, #> e.g. you will likely need to zoom in to popup-query features #> #> to see the number of points/vertices of the layer use 'npts(x)' #> to see the threshold for the feature type use 'mapview:::getMaxFeatures(x)' #> to adjust the threshold use argument 'maxpoints'
感想
どうでしょうか。
正直、メタデータの紐づけは思ったよりかなり難しかったです。ここ2週間ほど、元データのカオスっぷりに発狂しそうになりながら作業をしていました。死闘の日々は以下のレポジトリで垣間見れるので興味ある方はどうぞ。
実はまだきれいにデータを取り出しきれないんですけど、もう精神的に限界なのでこのへんで終わりにすることにしました。。e-Stat APIがつらいと思ってたんですけど、あれは恵まれてたんだなあとしみじみしてしまいました。
まとめ
データをスクレイピングするとき、あなたの精神もまたデータにスクレイプされているのだ。
