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

メッシュデータについて

GISデータでメッシュ単位でデータを公開している場合があるが、 メッシュデータにはいくつかの決まり事があるらしい。

詳しくは総務省統計局パスコのHPArcGIS Tipsがわかりやすい。

代表的なものとして第1次~第3次メッシュがあり、以下のようになっている。

名前 メッシュ単位 メッシュコード
第一次メッシュ 約80km四方 4桁
第二次メッシュ 約10km四方 6桁
第三次メッシュ 約1km四方 8桁

第一次メッシュは1度ごとの経線と2/3度ごとの緯線によって全国を分割して作られたものだそうである。それらを8等分したのが第2次メッシュ、さらに10等分して作られたのが第3次メッシュだそうである。また、メッシュ単位に とついているとおり、地域(緯度経度)によって多少メッシュの大きさが異なるそうである。

メッシュコードは緯度経度をもとに算出されているそうであるが、 ここでは割愛する。

市町村別メッシュコードは総務省統計局で公開している。

Rから国土数値情報をWebAPIを利用するkokudosuuchiパッケージが公開されている。

これらを使えば容易にデータが取得できそうである。 ただ、市町村別メッシュコードの更新年度が2010年なので少し注意が必要かもしれない。

mapviewパッケージによるインタラクティブな作図

leaflet::leaflet() で基本図の設定を行ったうえで、 mapview::garnishMap()関数を用いることにより様々な図を 重ねていくことができる。

最後にaddLayersControl()により、レイヤーの構成を調整することができる。

library(mapview)
library(sf)

nc <- st_read(system.file("shape/nc.shp", package = "sf"))
nc1 <- nc[1:50, ]
nc2 <- nc[51:100, ]
m_base <- 
  leaflet((options =leafletOptions(crs = leafletCRS(code="EPSG:4267"))) %>% 
  addTiles(group = "base") %>% 
  addProviderTiles(provider = "Stamen.TonerHybrid", group = "base_steman")

m <- garnishMap(m_base, addPolygons, data = as(nc2, "Spatial"), 
                weight = 2, group = "nc2")
m <- garnishMap(m, addPolygons, data = as(nc1, "Spatial"), 
                weight = 2, group = "nc2")

m <- addLayersControl(m, 
                      baseGroups = c("base", "base_steman"), 
                      overlayGroups = c("nc1", "nc2"))
m

f:id:jerrarrdan:20170423084755p:plain

ポリゴン内のラインデータの抽出(切取り)

何回も調べたりしているため、メモ。

ポリゴン内のライン長さを出したいときがある。 その際はラインをポリゴンでインターセクトする必要がある。

関数はrgeos::gIntersection()を利用する。

sp::over()関数で行うと意図と異なる結果となるため注意。
(適切な方法で行えば実行できるだろうが。。。)

library(stplanr)
library(magrittr)
data(routes_fast)
data(zones)
  • よい例
plot(zones[3,], col = "yellow")
plot(routes_fast, col ="grey90", add = T, lwd = 5)

rgeos::gIntersection(spgeom1 = routes_fast,
                     spgeom2 = zones[3,]) %>%
  plot(add = T, col = "blue")

f:id:jerrarrdan:20170418235229p:plain

  • ダメな例
plot(zones[3,], col = "yellow")
plot(routes_fast, col ="grey90", add = T, lwd = 5)

tmp <-
  sp::over(x = routes_fast, y = zones[3, ])

plot(routes_fast[!is.na(tmp$geo_code), ], add = T, col = "red",
     lwd = 1)

f:id:jerrarrdan:20170418235238p:plain

sf classデータの操作及びプロット

sfパッケージ(クラス)は空間データを簡単に操作することができ非常に便利である。

spパッケージ(クラス)はそれらを扱う多様なパッケージが存在するが、 自分自身で多少の編集を行うのは少し面倒であった。

Simple Features for Rだけあってデータ操作は簡単に行うことができる。 dplyr関係のものもサポートされており、 ggplot2パッケージではこれらのクラスを扱うgeom_sf()関数も開発されており、 今後も開発が活発に行われることが想像される。

簡単な例題を残しておく

  • mutate
nc <- sf::st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
nc %>% 
  mutate(A2 = AREA/ max(AREA)) %>%
  ggplot(data = .)+
  geom_sf(aes(fill = AREA))

f:id:jerrarrdan:20170417214211p:plain

  • filter
nc %>% 
  filter(AREA > 0.1) %>%
  ggplot(data = .) +
  geom_sf(aes(fill = AREA))

f:id:jerrarrdan:20170417214155p:plain

  • color gradient
nc  %>%
  ggplot()+
  geom_sf(aes(fill = AREA))+
  theme_bw() +
  scale_fill_distiller(palette = "YlOrRd", direction = 1)

f:id:jerrarrdan:20170417214136p:plain

  • annotate
ggplot(nc) +
  geom_sf() +
  annotate("point", x = -80, y = 35, colour = "red", size = 4)

ggplot(nc) +
  geom_sf() +
  annotate("point", x = -80, y = 35, colour = "red", size = 4) +
  geom_text(aes(x= -80, y = 35, label = "annotate"), hjust = "inward", vjust = "inward") 

f:id:jerrarrdan:20170417213756p:plain f:id:jerrarrdan:20170417214050p:plain

ggpairs()できれいなペアプロット

pairs()をggplot2風に書けるらしい。

library(tidyverse)
library(GGally)

iris %>% 
  GGally::ggpairs(data = ., aes(color = Species))

f:id:jerrarrdan:20170414232508p:plain

美しいペアプロット図を簡単に作る

rmarkdwonでの参考文献

ヘッダーは下記のようにする。

---
output:
  html_document:
  bibliography: bibdb.bib
  csl: japan-journal-of-industrial-and-applied-mathematics.csl
---

文章中の参考としたところには[@****]と表示してやればよい。

参考文献のスタイルを指定するcslファイルはZetoro style repositoryから自分に合うものを探す。

そして最終行を下記のようにする。(書かなくても表示はされる)

# References

参考したサイト

sfクラスからSpatial*****DataFrameクラスへの変換

Rにおけるジオメトリデータの操作は今後sfパッケージが主流になるらしい。 http://notchained.hatenablog.com/entry/2017/01/06/213333

開発が活発なためせっかく整備したコードが今後も使えるという保証はないため, 利用に少しためらいもあるが,最近少し使っている。

sfクラスからSpatial***DataFrameクラスに変換する際のメモである。

通常は以下のようにas(., "Spatial")で従来のSpatial****DataFrameに変換することができる。

    > nc <- st_read(system.file("shape/nc.shp", package="sf"), quiet = TRUE)

    > > class(nc)
    [1] "sf"         "data.frame"

    > methods::as(nc, "Spatial") %>% class
    [1] "SpatialPolygonsDataFrame"
    attr(,"package")
    [1] "sp"

しかし,tbl_dfクラスも持つ場合にはエラーが生じてしまう。 その場合は以下のようにして,クラスをsf, data.frameにすると解決できる。

    nc %>%
    `class<-`(c("sf", "data.frame")) %>%
    methods::as(., "Spatial")