Rでインタラクティブな地図表示はmapviewで十分

たまたまmapviewパッケージのヘルプ(example)をみていたら、 mapview()関数を + でつないでいく方法があることを知った。

もう単なる地図表示だけとりあえず表示してインタラクティブに確認したい場合は これを使えば十分そう。

ただし、データ容量の大きい場合には注意がいりそうなコメントをどこかで見たので、 その点は注意が必要そうである。

library(mapview)
library(kokudosuuchi)
library(rgdal)
library(sf)
library(ggplot2)
library(dplyr)
library(sp)

rgdal::setCPLConfigOption("SHAPE_ENCODING", "")

# データの取得

get_url <- kokudosuuchi::getKSJURL(identifier="A38", prefCode = 18)
## 医療圏
A38 <- kokudosuuchi::getKSJData(zip_url = get_url$zipFileUrl)
A38 <- A38[[1]]

## 医療機関
get_url <- kokudosuuchi::getKSJURL(identifier="P04", prefCode = 18)
P04 <- kokudosuuchi::getKSJData(zip_url = get_url$zipFileUrl)
P04 <- P04[[1]]

plot(A38, axes=T)
plot(P04, add=T, col="green")

A38 <- A38 %>% st_as_sf() %>% st_set_crs(4612) %>% st_transform(2448)
P04 <- P04 %>% st_as_sf() %>% st_set_crs(4612) %>% st_transform(2448)

m <- mapview(A38, zcol="A38a_004") + mapview(P04, zcol="診療科目1")

原因不明だが、このmapview()関数の順序を逆にするとうまく表示されなかった。

mapshot(x = m, file = "fukui_iryou.png")

ついでに、医療圏域別に診療科目の開設？状況について表にしてみた。 人口別にみる必要があると覆うがかなり偏りのある分布になった。 徒歩圏内に医療がない地域も多そうで、医療の問題は大変そう。

歯科が多いのはよく聞くが、内科や小児科ってかなり多いんですね。 安心して子供を産める地域といえるのかもしれませんが、 小児科の数は少子化社会に対して妥当な数なんでしょうかね。

診療科目1に複数の科目が入力されているため切り分けて集計する。

st_intersection(x = P04, y = A38) %>%
  as.data.frame() %>%
  select(診療科目1, A38a_004) %>%
  tidyr::separate(col = 診療科目1, into = paste0("診療科目", 1:23), sep = "　") %>%
  mutate(rn = row.names(.)) %>% 
  gather(key = 科目, value = prof, starts_with("診療科目")) %>%
  xtabs(~A38a_004+prof, data = .) %>%
  as.data.frame() %>%
  spread(A38a_004, Freq) %>%
  mutate(計 = `奥越` + `丹南` + `福井・坂井` + `嶺南`) %>%
  arrange(desc(`計`))

prof	奥越	丹南	福井・坂井	嶺南	計
内科	25	93	191	62	371
歯科	21	64	158	46	289
小児科	15	56	113	34	218
消化器科	13	35	83	21	152
外科	10	39	64	23	136
整形外科	9	33	74	19	135
リハビリテーション科	8	27	70	15	120
循環器科	8	27	72	11	118
皮膚科	7	25	45	17	94
呼吸器科	6	18	50	7	81
眼科	4	14	35	9	62
耳鼻咽喉科	4	11	31	11	57
放射線科	3	21	25	8	57
小児歯科	0	12	31	6	49
肛門科	5	12	26	5	48
リウマチ科	2	7	34	4	47
泌尿器科	2	13	22	7	44
アレルギー科	4	2	30	5	41
胃腸科	2	12	19	4	37
精神科	1	4	22	9	36
矯正歯科	0	7	25	2	34
脳神経外科	1	9	18	4	32
産婦人科	2	6	16	6	30
神経内科	2	5	19	4	30
歯科口腔外科	0	4	20	5	29
麻酔科	2	7	15	5	29
形成外科	1	6	15	2	24
婦人科	1	8	9	5	23
心療内科	0	3	14	3	20
神経科	0	3	8	5	16
気管食道科	1	0	11	1	13
心臓血管外科	0	3	7	3	13
産科	0	3	6	1	10
呼吸器外科	0	1	4	3	8
性病科	0	4	1	0	5
美容外科	0	1	2	0	3
小児外科	0	0	2	0	2
皮膚泌尿器科	0	1	0	0	1

歯科はコンビニより多いということをどこかのTV番組 できいたので、多いらしいということは知っていたけど、それを上回って内科が多い。 ちなみに福井県の場合、コンビニ数はサイト確認時で353店舗なので、内科の病院のみがコンビニ店舗数を上回っている。

ポリゴン内にテキスト（ラベル）を表示する

ポリゴン内に各ポリゴンの属性を表示したい。 何回か調べているので、メモとして残す。

# using packages
library(ggplot2)
library(kokudosuuchi)
library(sf)
library(sp)
options(stringsAsFactors = F)

# preparation
url <- kokudosuuchi::getKSJURL("N03", prefCode = c(27))
url <- url %>% arrange(desc(year)) %>% .[1, ]
N03 <- kokudosuuchi::getKSJData(url$zipFileUrl)
N03 <- N03[[1]]
proj4string(N03) <- CRS("+init=epsg:4612")

# sp class
plot(N03, axes=T)
text(getSpPPolygonsLabptSlots(N03), labels = as.character(N03$市区町村名), cex=0.4)

# sf class
N03_sf <- st_as_sf(N03)
N03_sf %>% 
  st_geometry() %>%
  plot()
text(st_coordinates(N03_sf %>% st_centroid), labels=N03$市区町村名, cex = 0.3)

# using ggplot2
N03_sf %>%
  mutate(lon = st_centroid(.) %>% st_coordinates() %>% .[, 1],
         lat = st_centroid(.) %>% st_coordinates() %>% .[, 2]) %>%
  ggplot() +
  geom_sf(aes(fill = 市区町村名)) +
   ggrepel::geom_text_repel(aes(x=lon, y=lat, label=市区町村名), size=2)

こちらはstackoverflowにあった回答例からのコピペ。 ggplot2パッケージでもっと簡単にテキストを作図したいが、わからず。

ここ一年のTOPIX30の株価（終値）の推移

ここ一年のTOPIX30の株価（終値）を可視化してみた。

便利なパッケージがあるらしいが、うまく機能しなかった。 yahooではスクレイピングを禁止しているらしい、 また、させないように仕様を頻繁に？変更しているらしいので、 そのほかでできそうな記事を参考にやってみた。

これだけで各株のトレンドが一目瞭然なので便利そう。

library(tidyverse)
library(lubridate)
options(stringsAsFactors = F)

# 証券コードと名前の変換用
code_list <- 
  rio::import("http://www.jpx.co.jp/markets/statistics-equities/misc/tvdivq0000001vg2-att/data_j.xls")

topix30 <- code_list %>% filter(規模区分 == "TOPIX Core30")
code <- topix30$コード

stock <- NULL
for (i in seq_len(length(code))){
  url <- paste("http://k-db.com/stocks/", 
               code[i], 
               "-T?1d/2017?&download=csv", 
               sep = "")
  stock[[i]] <- read.table(url, sep = ",", skip = 2)
  colnames(stock[[i]]) <- 
    c("date", "start_price", "max_price", "min_price",
      "end_price","volume","trading_value")
  # dateをでdate要素に変更
  stock[[i]] <- stock[[i]] %>% mutate_at(.vars = "date", ymd)
  Sys.sleep(5) 
}
names(stock) <- code

data_list <- NULL
for (i in seq_len(length(stock))){
  data_list <- rbind(data_list, 
                     stock[[i]] %>% mutate(code = as.character(code[i])))
}

p <- 
  data_list %>%
  mutate(code = plyr::mapvalues(x = code, from = topix30$コード, to = topix30$銘柄名)) %>% 
  ggplot() +
  geom_line(aes(x = date, y = end_price)) +
  facet_wrap(~code, scales = "free_y") 

p

点は月曜日である。 JR系の株価がうなぎのぼり。

データの取得に関してはほとんどこのサイト のパクリです。

自作関数の高速化

efficient R programingを読んでいて今後使えそうなことに関するメモ。

3.7 では自作関数の高速化について述べられている。

自作関数にcompiler::cmpfun()を適用すればよいだけらしい。

library(compiler)

mean_r = function(x) {
  m = 0
  n = length(x)
  for(i in seq_len(n))
    m = m + x[i] / n
  m
}

cmp_mean_r = cmpfun(mean_r)

平均を求める関数についての比較結果はこのようになる。 結構差が出ているので、複雑な処理をする際には有効に使えそうである。

x = rnorm(1000)
microbenchmark(times = 10, unit = "ms", # milliseconds
          mean_r(x), cmp_mean_r(x), mean(x))
#> Unit: milliseconds
#>           expr   min    lq  mean median    uq  max neval cld
#>      mean_r(x) 0.358 0.361 0.370  0.363 0.367 0.43    10   c
#>  cmp_mean_r(x) 0.050 0.051 0.052  0.051 0.051 0.07    10  b 
#>        mean(x) 0.005 0.005 0.008  0.007 0.008 0.03    10 a  

• 2017/7/23追記 R言語逆引きハンドブックにバイトコンパイルによる高速化として、 しっかり書いてあったので、詳細はそちらを確認。

tmapとggplot2の出力結果の比較

空間データのプロット

sfやspクラスオブジェクトのデータをプロットする方法について tmapとggplot2パッケージによる出力を比較してみる。

library(tmap)
library(sf)
library(mapview)
library(ggplot2)

spクラスオブジェクトとsfクラスオブジェクトの2つを準備する。

nc <- st_read(system.file("shape/nc.shp", package="sf"), quiet = T)
data("atlStorms2005")

class(nc)

[1] "sf"         "data.frame"

class(atlStorms2005)

[1] "SpatialLinesDataFrame"
attr(,"package")
[1] "sp"

tmapを利用して書くと

tm_shape(nc) + tm_polygons(col = "AREA")

tm_shape(atlStorms2005) + tm_lines(col = "MaxWind") +
  tm_layout(frame = F)

ggplot2を利用して書くと

ggplot() +
  geom_sf(data = nc, aes(fill = AREA))

ggplot() +
  geom_sf(data = st_as_sf(atlStorms2005), aes(color = MaxWind))

tmapパッケージではsp, sfオブジェクトともに区別なく書くことができるが、 ggplot2ではsfオブジェクトのみ対応しており、変換が必要になる。

ggplot2パッケージのパレットが見づらいと感じるので、 デフォルトのパレットを変更されることを期待している。

行の結合にはrbind.fill()が便利

各データフレームを行で結合したい場合がある。 列で結合する場合はleft_join()などのjoin関数がある。

データがない部分はNAで保管される。

    library(plyr)
    rbind.fill(mtcars[c("mpg", "wt")], mtcars[c("wt", "cyl")])

    mpg    wt cyl
    1 21.0 2.620  NA
    2 21.0 2.875  NA
    3 22.8 2.320  NA
    4 21.4 3.215  NA
    5 18.7 3.440  NA
    6 18.1 3.460  NA
    .................
    59  NA 2.140   4
    60  NA 1.513   4
    61  NA 3.170   8
    62  NA 2.770   6
    63  NA 3.570   8
    64  NA 2.780   4

• 2017/7/23追記 dplyrパッケージのbind_rows()、rbind_all()関数でも同様のことができるらしいので、 わざわざplyrパッケージをロードして使用する必要もないらしい。 ただし、bind_rows()はオブジェクトクラスによってはうまく機能しないらしいので注意。 現時点ではsfオブジェクトに対してはエラーが発生する。

leaflet、mapviewパッケージを用いた可視化のコードサンプル

いつも調べながら書いているので、 よく利用するオプションを書いておく

# 必要とするパッケージ
library(mapview)
library(kokudosuuchi)
library(tidyverse)
library(sp)

# データの読込(岐阜県の行政区域)
getURL <- 
  kokudosuuchi::getKSJURL("N03", prefCode=21) %>% 
  filter(year == 2017) %>%
  .[,"zipFileUrl"]

adm_area <- kokudosuuchi::getKSJData(as.character(getURL))
adm_area <- adm_area[[1]]
proj4string(adm_area) <- CRS("+init=epsg:4612")

# mapview,leafletパッケージを用いて可視化
pal <- colorFactor(palette = "YlOrRd", domain = adm_area$市区町村名)
leaflet() %>%
  addProviderTiles("OpenStreetMap.BlackAndWhite") %>%
  addPolygons(data = adm_area, weight = 1, group = "adm area", 
              fillColor = ~ pal(adm_area$市区町村名),
              fillOpacity = 0.5,
              highlightOptions = highlightOptions(color = "green", weight = 7),
              popup = popupTable(adm_area@data)) %>%
  addLayersControl(overlayGroups = c("adm area"), position = "bottomright") %>%
  # mapview::mapshot(url = "test.html") # 静的な図として保存する場合はコメントアウト