読者です 読者をやめる 読者になる 読者になる

まだ厨二病

RとPyhtonを使ったデータ分析・統計解析らへんの話題をしていくだけ

⚡️関数の実行結果を知らせる関数(パッケージ)を作っています

source()でガーッとコードを実行させる時やforeach::foreach()などでループ処理をする時、モデリングのような処理に時間のかかる関数を実行すると、終わったかな、まだやっているかなとチラチラとコンソールを確認するのが結構面倒臭い。せっかちなので頻繁に見てしまう。そんな私のための、関数の実行が終了したら通知する関数を書きました。パッケージとしてGitHubに公開しているので、後述の方法でインストールすれば利用可能です。

まずはこちらをご覧ください。

f:id:u_ribo:20170114083229g:plain

このgifは下記のコードの実行結果を記録したものです。

library(roracle)
notice_status(head(iris))

head(iris) %>% notice_status(msg = "実行完了")

notice_status(warning("警告の際には出力が変更されます"))
notice_status(print(エラーとなるコードを実行するとそれを通知します))


notice_overtime(1 + 1)
notice_overtime(for(i in 1:1e7) sqrt(1:10)) 

インストール

{roracle}というパッケージを利用することで上のような通知が可能になります。下記のコードを実行すればGitHubからインストールされます(さっきGitHubにあげたのでGepuro Task Viewsに登録されていないはず)。

github.com

install.packages("devtools")
install_github("uribo/roracle")

通知のためにnotiというツールをインストールしておく必要があります。これが{roracle}の味噌となります(つまり私自身は大したことをしていない!!)

github.com

notiはOSを問わずに利用できます。

使い方

{roracle}の関数を紹介します。といっても、現在は2つの関数しかありません。コードの実行結果を通知するnotice_status()と処理時間の長いコードの終了を知らせるnotice_overtime()です。

それぞれ詳しく見ていきます。

notice_status()では、第一引数に与えたコードが実行されると通知をします。%>%によるパイプ処理をしても良いです。msgという引数を持っていて、これを変更すると通知のメッセージを変えられます。またそのうち音声読み上げさせたいと思っています。

また先の例のように、コード実行結果にWarningやMessage、Errorがあった場合にはその旨を通知します。

notice_status(head(iris))

notice_overtime()は実行時間が規定以上になったものに対して、実行完了時に通知を行います。s処理時間の長い関数の実行に追加しておくと、終了したことを確認できて捗ります。

これで一つの画面を見続けている必要がなくなりました。

Enjoy!

この関数を書いた背景として以前の記事があります。

uribo.hatenablog.com

uribo.hatenablog.com

ただこれだとMacでしか動かなかったり、パッケージ化されていないので使いにくい、ということがありました。個人的には満足できるものができたと思います。CRANにあげる気はありません。関数を使わず、もっと自然に通知させたいなと思っていますが、それはまた今度。

また、同様のことをしようとしている人がいます。metacranおじさん(Mango Solutionsの人)もなかなかの変態です。いつもありがとう!!

Rで植生図を描画する

植生図というものがあります。利用目的に応じた特定の項目(主題)を表現した主題図の一つです。、主題図には他に人口の内訳とか土地利用図などがあります。植生図は植生、すなわち、ある地域に生育する植物集団の特徴を地図上に表現したものと思ってもらって構いません。

さて、日本全国の植生の現状を整備、保護・活用することを目的とした調査が環境省によって行われています。これは「自然環境保全基礎調査」の中の「植生調査」と呼ばれる調査なのですが、過去数回の調査の成果物として植生図GISデータおよび植生図画像が提供されており、現在も最新の調査結果である平成26年度・27年度のものが追加・更新されています。

例えば2.5万地形図 図葉名「箭田」ではこちらのページのリンクにあるような画像ファイルです。前置きが長くなりましたが今回は、すでに公開されている植生図について、Rで作ってみようという話です。

用意するもの

方法

ゆたにーこと、id:yutannihilationによるとRでGISをやるときにはsfパッケージになるようなのですが、古い人間なので{maptools}Shapefileを読み込みます。勉強不足による敗北感があります{sf}を使った地理空間データの読み込みではsfというオブジェクトとして扱われますが、{maptools}では従来の地理空間クラスの一種であるSpatialPolygonsDataFrameとして扱われます。(余談。{sf}には未来を感じます

# devtools::install_github("tidyverse/tidyverse")
library(tidyverse)
library(forcats)
library(magrittr)

# 513375 箭田のShapefileファイルを読み込む
map <- maptools::readShapeSpatial("~/Downloads/vg67_33/shp513375/p513375.shp", 
                                  proj4string = sp::CRS("+proj=utm +zone=53 +ellps=WGS84 +units=km"))
# dataの変数名を確認
map@data %>% names()
##  [1] "MESH2_C"   "HANREI_C"  "SURV_YEAR" "ORG_NO"    "ZUKAKU_NO"
##  [6] "SHOKU_C"   "SHOKU_N"   "DAI_C"     "DAI_N"     "CHU_C"    
## [11] "CHU_N"     "SAI_C"     "SAI_N"     "HANREI_N"
# 文字列が文字化けしているので修正
map@data$HANREI_N[1]
## [1] \x90\x85\x93c\x8eG\x91\x90\x8cQ\x97\x8e
## 30 Levels: \x83\x84\x83i\x83M\x8d\x82\x96،Q\x97\x8e\x81i\x82u\x82h\x81j ...
map@data %<>% mutate_if(is.character, stringr::str_conv, encoding = "SJIS")
map@data$HANREI_N[1]
## [1] "水田雑草群落"

ここで示したHANREI_Nという変数が、植生調査の結果として得られた統一凡例名となります。箭田データでは、合計で30の凡例が使用されています。

map@data$HANREI_N %>% unique() %>% length()
## [1] 30

凡例の調整

植生図では、多様な群落を色や塗り分けの種類を分けることが大事です。自然環境保全基礎調査の植生図でも、凡例に応じた塗り分けが実施されています。ですが、元のファイルに凡例と対応した色の指定がないので、ここはちまちまと画像ファイルを見ながらそこで使われている色を凡例番号として利用されるORG_NOという列と対応させる処理をとります。

map@data %<>% mutate(
                     color = case_when(
                       .$ORG_NO == "12" ~ "#DBC9AAFF",
                       .$ORG_NO == "16" ~ "#FAE6A6FF", 
                       .$ORG_NO == "17" ~ "#FEE085FF", 
                       .$ORG_NO == "18" ~ "#FED9CDFF", 
                       .$ORG_NO == "2" ~ "#B4E983FF", 
                       .$ORG_NO == "20" ~ "#BFE5DEFF", 
                       .$ORG_NO == "21" ~ "#98F7EDFF", 
                       .$ORG_NO == "22" ~ "#C6C5F3FF", 
                       .$ORG_NO == "23" ~ "#F8D9DCFF",
                       .$ORG_NO == "3" ~ "#95DE82FF",
                       .$ORG_NO == "4" ~ "#C6B8B8FF",
                       .$ORG_NO == "5" ~ "#C28182FF",
                       .$ORG_NO == "6" ~ "#A67F81FF",
                       .$ORG_NO == "7" ~ "#80C0E1FF",
                       .$ORG_NO == "8" ~ "#9981BEFF", 
                       .$ORG_NO == "9" ~ "#FFF585FF", 
                       .$ORG_NO == "a" ~ "#FFF2D9FF",
                       .$ORG_NO == "b" ~  "#C7ECFEFF", 
                       .$ORG_NO == "c" ~ "#F3D8CBFF",
                       .$ORG_NO == "d" ~  "#B3F3F2FF",
                       .$ORG_NO == "e" ~  "#F3D8CCFF",
                       .$ORG_NO == "f" ~   "#F3FDE3FF", 
                       .$ORG_NO == "g" ~ "#FEFD87FF",
                       .$ORG_NO == "h" ~   "#F9FD87FF", 
                       .$ORG_NO == "i" ~ "#F2F2F2FF", 
                       .$ORG_NO == "k" ~ "#CCD3CCFF", 
                       .$ORG_NO == "L" ~  "#BFBFBFFF",
                       .$ORG_NO == "m" ~ "#BAACA9FF", 
                       .$ORG_NO == "r" ~ "#FED185FF", 
                       .$ORG_NO == "w" ~ "#E6FFFCFF"))

次に凡例を自然環境保全基礎調査のもので使われている順番に並び替えます。一度、凡例番号と統一凡例名からなる文字列を要因データとし、forcats::fct_relevel()により水準の並びを入れ替えます。

map@data %<>% mutate(hanrei_new = fct_relevel(
  as.factor(paste(ORG_NO, HANREI_N)),
  c("3 カナメモチ-コジイ群集", "22 ヤナギ高木群落(VI)", "4 シイ・カシ二次林", "2 アベマキ-コナラ群集", 
    "5 アカマツ群落(VII)", "6 ネズ-アカマツ群落", "17 低木群落", "23 クズ群落", 
    "9 ススキ群団(VII)", "18 伐採跡地群落(VII)", "8 ヌマガヤオーダー", "7 ヨシクラス", 
    "20 ツルヨシ群集", "21 オギ群集","12 スギ・ヒノキ・サワラ植林", "16 竹林",
    "h ゴルフ場・芝地", "g 牧草地", "f 路傍・空地雑草群落", "c 放棄畑雑草群落", "e 果樹園", 
    "a 畑雑草群落", "b 水田雑草群落", "d 放棄水田雑草群落", "k 市街地", 
    "i 緑の多い住宅地", "L 工場地帯", "m 造成地", "w 開放水域", "r 自然裸地")))

これで準備は完了です。

主題図を描く

肝心の地図描画作業です。ここでは{tmap}というパッケージを使って楽に主題図を描画させることができます。{tmap}はRにおける定番グラフ描画ライブラリの{ggplot2}をベースとしており、{ggplot2}のようにレイヤーを重ねていくことで目的の図を作成していきます。また日本語のラベルを使うために{extrafont}も読み込んでおきます。

library(tmap)
library(extrafont)
# フォントの読み込み
loadfonts(quiet = TRUE)

tm_shape()tm_fill()を使って地図とデータを描画させたものをベースとし、ここに要素を重ねていきます。出来る限り元の植生図を再現したものを作っていきましょう。

colorsets <- map@data %>% arrange(hanrei_new) %>% use_series(color) %>% unique()

vg.map <- tm_shape(map) +
  tm_fill("hanrei_new", title = "植生図凡例", palette = colorsets, alpha = 0.6)

凡例の位置や、縮尺、方位記号、タイトルや各種の配置を調整し、完成です。

vg.map +
  tm_text("ORG_NO", size = 0.35) +
  tm_borders(lwd = 0.6, col = "black", alpha = 0.5) +
  # 縮尺の追加
  tm_scale_bar(position = c("right", "bottom")) +
  # 方位記号の追加
  tm_compass(type = "4star", color.dark = "red", show.labels = 2, size = 2, position = c("right", "top")) +
  tm_layout("513375 箭田", 
            title.position = c("left", "top"), 
            title.size = 1,
            frame = FALSE, 
            fontfamily = "IPAexGothic",
            legend.outside = TRUE,
            inner.margins = c(0.10, 0.06, 0.10, 0.12)) +
  tm_credits("1/25,000植生図「箭田」GISデータ(環境省生物多様性センター)を使用し、瓜生真也が作成・加工したものである。\n(http://gis.biodic.go.jp/webgis/sc-002.html#webgis/513375)",
             size = 0.4,
             position = c("left", "bottom")) +
  tm_legend(bg.color = "white", 
            # スケールバーを枠外に
            # attr.outside = TRUE,
            # legend.outside.size = 0.2,
            legend.text.size = 0.6,
            frame = FALSE, legend.outside = TRUE)

f:id:u_ribo:20170107063100p:plain

いい感じです。凡例と色の紐付けだけが面倒ですが、簡単に植生図を描画できました。

おまけ: leafletでも

{leaflet}で植生図を描画させると、拡大して確認できるしタイルを重ねられるので便利ですね。やっていきましょう。せっかくなので国土地理院提供の地理院タイルを背景に描画させます([参考](http://notchained.hatenablog.com/entry/2015/10/26/231508)。

library(leaflet)
atr <- "<a href='http://maps.gsi.go.jp/development/ichiran.html' target='_blank'>地理院タイル</a>"

leaflet(map) %>%
  addTiles("http://cyberjapandata.gsi.go.jp/xyz/pale/{z}/{x}/{y}.png",
           attribution = atr) %>%
  addPolygons(
    stroke = FALSE, 
    fillOpacity = 0.8, 
    smoothFactor = 0.5,
    color = ~color,
    label = ~HANREI_N,
    # labelOptionsはleaflet開発版(1.0.2.9008)の機能です
    labelOptions = labelOptions(opacity = 1, 
                   noHide = TRUE, 
                   textOnly = FALSE,
                   direction = "auto", offset = c(20,-15))
  )

f:id:u_ribo:20170107064001g:plain

Enjoy!

出典

このページで表示した植生図は、1/25,000植生図「箭田」GISデータ(環境省生物多様性センター)を使用し、瓜生真也が作成・加工したものである。(http://gis.biodic.go.jp/webgis/sc-002.html#webgis/513375)

⭐️Rを使ったモデル構築の最善策を求めて: {dplyr} + {tidyr} + {broom} + {purrr}を使ったアプローチ

RStudioのチーフサイエンティスト、Hadley Wickham(ハドリー)が2月に行った講演のビデオがYouTubeに上がっていたので観た。

"Making Data Analysis Easier"というタイトルでの発表(スライドでは"Managing many models"になっているけど)で、ハドリー自身が考えている、データサイエンスに必要な可視化やモデリングを効率的に行うための手法について、彼の開発してきたパッケージを中心に説明している。

www.youtube.com

分かりやすく、具体例を交えた内容なので、是非YouTubeの動画を観てもらうのが良いと思うが、自分の頭を整理するためにもここでまとめておく。なお、発表スライドはクリエイティブ・コモンズライセンス3.0のもと、表示・非営利のラインセンスで再利用可能となっている。

Hadley Wickham (Chief Scientist, RStudio) 2016. Managing many models. http://wombat2016.org/abstracts/hadley.html

  • ✨ データサイエンスにおいて重要な7つの事柄
  • 🔰 データの階層構造を理解する
    • 現実的な問題: gapminderデータを例に
    • ある基準でデータを内包する
    • {broom}パッケージでRによるモデルの結果をよしなに
    • {dplyr}と{purrr}および{broom}による効率的なモデリング
  • ⭐ まとめ
  • 💻 実行環境

✨ データサイエンスにおいて重要な7つの事柄

ハドリーは自身の発表でしばしば次のような図を提示する。

f:id:u_ribo:20160326224249p:plain

続きを読む